Host-Parasite Meta-analysis
Amanda Vicente, Beatriz Willink, Kacy Nowak, Dave Civitello
2023-06-22
Load packages
x <-
c(
"metafor", "kableExtra", "Matrix", "dplyr", "wesanderson", "weightr", "ggplot2", "MuMIn", "RColorBrewer", "tidyverse", "patchwork", "devtools", "R.rsp", "orchaRd")
lapply(x, function(y) {
# check if installed, if not install
if (!y %in% installed.packages()[, "Package"])
install.packages(y)
# load package
try(require(y, character.only = T), silent = T)
})
#remotes::install_github("itchyshin/orchard_plot", subdir = "orchaRd", force = TRUE, build_vignettes = TRUE)Read in data
dat <- read.csv("../data/MA_dat-21-06-2023.csv", header = T, sep = ",")
# subset data that have a "variation type". We'll parse based on variation type below
#dat <- dat[-which(is.na(dat$Variation.Type) == TRUE),]
dat <- subset(dat, select = -c(Data_source,Comment))
colnames(dat)## [1] "ID" "Study"
## [3] "Experiment" "Entered.by"
## [5] "Host" "Host.type"
## [7] "Parasite" "Parasite.type"
## [9] "Study.type" "X_gradient"
## [11] "Gradient.category" "Trait"
## [13] "Trait.type" "Trait.category"
## [15] "Gradient.as.predicted" "Statistics"
## [17] "Level_C" "Level_X"
## [19] "N_C" "N_X"
## [21] "Mean_C" "Variation_C"
## [23] "Lwr_C" "Upr_C"
## [25] "Mean_X" "Variation_X"
## [27] "Lwr_X" "Upr_X"
## [29] "Success_C" "Success_X"
## [31] "Variation.Type" "OR"
## [33] "OR_SE" "OR_N"
## [35] "F_X" "F_df"
## [37] "Chi2" "Chi2_df"
## [39] "Slope" "R"
## [41] "R_N" "R_df"
## [43] "Slope_SE" "Z_X"
## [45] "Z_N" "Z_df"
## [47] "Z_SE" "Environment"
## [49] "Environmental_stage" "Vector_transmitted"
## [51] "Final_host" "Require_intermediate_host"
## [53] "DOI"Data cleaning
Gradient categories
# Toxin is a type of pollution
dat$Gradient.category <- gsub(pattern = "Toxin", replacement = "Pollution", x = dat$Gradient.category)
# exclude ambiguous/equivocal categories
dat <- dat[which(dat$Gradient.category != "Ecological"),]
# Check gradient categories
levels(as.factor(dat$Gradient.category))## [1] "Environment" "Pollution" "Resource"# Check pollution experiments are actually chemical pollution
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Pollution"]))## [1] "17-estradiol (E2)"
## [2] "Alexandrium tamarense toxin"
## [3] "alpha-cypermethrin"
## [4] "Atrazine"
## [5] "Atrazine Herbicide"
## [6] "Cadmiun"
## [7] "Carbaryl"
## [8] "Carbaryl Pesticide "
## [9] "Copper contamination"
## [10] "Deltamethrin"
## [11] "diatomaceous earth"
## [12] "Diatomaceous earth"
## [13] "Diatomaceous Earth (DE)"
## [14] "Fipronil"
## [15] "Fipronil Insecticide"
## [16] "Flupyradifurone"
## [17] "Funguscides"
## [18] "Glyphosate"
## [19] "Glyphosate exposure"
## [20] "Imidacloprid"
## [21] "Insecticide"
## [22] "Lead Exposure"
## [23] "Malathion Insecticide Exposure "
## [24] "Naphthenic acids"
## [25] "Noctural light"
## [26] "Noise"
## [27] "Pesticide"
## [28] "pesticide lambda (k)-cyhalothrin"
## [29] "Pesticides (thiacloprid and tau-fluvalinate)"
## [30] "Thiacloprid"
## [31] "Thiacloprid Insecticide"
## [32] "Thiamethoxam"
## [33] "Wastewater treatment plant effluents"# Nocturnal light and noise should be "Environment"
dat$Gradient.category[dat$X_gradient == "Noctural light"] <- "Environment"
dat$Gradient.category[dat$X_gradient == "Noise"] <- "Environment"
# Check Environment experiments are not chemical pollution
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Environment"]))## [1] "Artificial light at night" "Diurnal Temperature Range "
## [3] "Drought" "Environmental enrichment"
## [5] "Humidity" "Humidity and Temperature "
## [7] "Hypoxia" "Light exposure"
## [9] "Noctural light" "Noise"
## [11] "Oxygen" "pH"
## [13] "Plant species diversity" "Salinity"
## [15] "Seasonality (Temperature) " "Temperature"
## [17] "tidal stress" "Water Temperature"# Check resource limitation experiments are actual resources
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Resource"]))## [1] "Food" "Food "
## [3] "Food availability" "Food level"
## [5] "Low food" "Nitrogen and phosphorus conc"
## [7] "Nutrional stress" "Pollen limitation"
## [9] "Protein (yeast) limitation" "Selenium supplementation"
## [11] "Starvation" "Starvation time"
## [13] "Sugar limitation"Transmission mode
Following R1 suggestion to discuss and/or include as moderator.
for (i in 1:nrow(dat)) {
if (is.na(dat$Environmental_stage[i]) == F &
is.na(dat$Vector_transmitted[i]) == F &
is.na(dat$Require_intermediate_host[i]) == F) {
if (dat$Environmental_stage[i] == "Present" |
dat$Vector_transmitted[i] == "Yes" |
dat$Require_intermediate_host[i] == "Yes") {
dat$Tansmission[i] <- "Indirect"
} else {
dat$Tansmission[i] <- "Direct"
}
} else {
dat$Tansmission[i] <- NA
}
}Median, min, max and n to mean and SD
non_par <- dat[which(dat$Variation.Type == "Median_Min_Max"),]
# Estimate mean for control and treatment from median, maximum and minimum
non_par$Estimated_Mean_C <- (non_par$Lwr_C + 2 * non_par$Mean_C + non_par$Upr_C) / 4
non_par$Estimated_Mean_X <- (non_par$Lwr_X + 2 * non_par$Mean_X + non_par$Upr_X) / 4
# Estimate SD from range and n
# First, read in Xi_N table from Wang et al. (2014) BMC Medical Research Methodology
Xi_table <- read.csv("../scripts/Xi_N_Table.csv", header = T, sep = ",")
# Estimate SD using the approximation for each N under 50
for (i in 1:length(non_par$ID)) {
non_par$Estimated_SD_C[i] <-
(non_par$Upr_C[i] - non_par$Lwr_C[i]) / Xi_table$Xi_N[which(Xi_table$N == non_par$N_C[i])]
non_par$Estimated_SD_X[i] <-
(non_par$Upr_X[i] - non_par$Lwr_X[i]) / Xi_table$Xi_N[which(Xi_table$N == non_par$N_X[i])]
}
# overwrite dat
dat$Mean_C[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_Mean_C
dat$Mean_X[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_Mean_X
dat$Variation_C[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_SD_C
dat$Variation_X[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_SD_X
dat$Variation.Type[which(dat$Variation.Type == "Median_Min_Max")] <- "SD"
# Sanity check
levels(as.factor(dat$Variation.Type))## [1] "" "CI" "Cloppler-Pearson CI"
## [4] "GLMM" "GM_SD" "IQ"
## [7] "logistic95CI" "Median_IQ" "OR"
## [10] "SD" "SE" "t"
## [13] "Wald CI" "Z" "Z_reg"Median and IQ range to mean
non_par_2 <- dat[which(dat$Variation.Type == "Median_IQ"),]
# Estimate mean for control and treatment from median and IQ range
non_par_2$Estimated_Mean_C <- (non_par_2$Lwr_C + non_par_2$Mean_C + non_par_2$Upr_C) / 3
non_par_2$Estimated_Mean_X <- (non_par_2$Lwr_X + non_par_2$Mean_X + non_par_2$Upr_X) / 3
# overwrite dat
dat$Mean_C[which(dat$Variation.Type == "Median_IQ")] <- non_par_2$Estimated_Mean_C
dat$Mean_X[which(dat$Variation.Type == "Median_IQ")] <- non_par_2$Estimated_Mean_X
dat$Variation.Type[which(dat$Variation.Type == "Median_IQ")] <- "IQ"
# Sanity check
levels(as.factor(dat$Variation.Type))## [1] "" "CI" "Cloppler-Pearson CI"
## [4] "GLMM" "GM_SD" "IQ"
## [7] "logistic95CI" "OR" "SD"
## [10] "SE" "t" "Wald CI"
## [13] "Z" "Z_reg"Subset effects with workable “variation types”
nrow(dat)## [1] 1194accepted_var <- c("SE", "SD", "CI", "Wald CI", "IQ", "OR")
dat <- dat[which(dat$Variation.Type %in% accepted_var),]
nrow(dat)## [1] 1144# store sample size in an object
N <- nrow(dat)Getting effect sizes
Calculate OR and approximate variance from contingency tables
for (i in 1:N) {
# no correction need if none of the categories are 0
if (dat$Variation.Type[i] == "OR" &
dat$Success_C[i] > 0 &
dat$Success_X[i] > 0 &
(dat$N_C[i] - dat$Success_C[i]) > 0 &
(dat$N_X[i] - dat$Success_X[i]) > 0) {
dat$OR[i] <-
(dat$Success_X[i] * (dat$N_C[i] - dat$Success_C[i])) / (dat$Success_C[i] * (dat$N_X[i] - dat$Success_X[i]))
dat$Log.OR[i] <- log(dat$OR[i])
# approximate variance
dat$Log.OR.v[i] <-
1 / dat$Success_X[i] + 1 / (dat$N_X[i] - dat$Success_X[i]) + 1 / dat$Success_C[i] + 1 / (dat$N_C[i] - dat$Success_C[i])
}
else {
# if at least one category is 0, we apply Yate's correction to avoid dividing by 0
if (dat$Variation.Type[i] == "OR" &
(dat$Success_C[i] == 0 |
dat$Success_X[i] == 0 |
(dat$N_C[i] - dat$Success_C[i]) == 0 |
(dat$N_X[i] - dat$Success_X[i]) == 0)) {
dat$OR[i] <-
((dat$Success_X[i] + 0.5) * (dat$N_C[i] - dat$Success_C[i] + 0.5)) / ((dat$Success_C[i] + 0.5) * (dat$N_X[i] - dat$Success_X[i] + 0.5))
dat$Log.OR[i] <- log(dat$OR[i])
# approximate variance
dat$Log.OR.v[i] <-
1 / (dat$Success_X[i] + 0.5) + 1 / (dat$N_X[i] - dat$Success_X[i] + 0.5) + 1 / (dat$Success_C[i] + 0.5) + 1 / (dat$N_C[i] - dat$Success_C[i] + 0.5)
}
else{
# if there is no OR data, make these columns NA
dat$OR[i] <- NA
dat$Log.OR[i] <- NA
dat$Log.OR.v[i] <- NA
}
}
}Calculate SD for comparisons with continuous normally distributed variables
# create SD variables
dat$SD_C <- vector(length = N)
dat$SD_X <- vector(length = N)
for (i in 1:N) {
# Calculate SD if SE is reported
if (dat$Variation.Type[i] == "SE") {
dat$SD_C[i] <- dat$Variation_C[i] * sqrt(dat$N_C[i])
dat$SD_X[i] <- dat$Variation_X[i] * sqrt(dat$N_X[i])
} else {
# calculate SD if a 95% confidence interval for a normal distribution is reported
# this also applies to the Wald confidence interval for proportions, as it is a normal approximation to the binomial
if (dat$Variation.Type[i] == "CI" |
dat$Variation.Type[i] == "Wald CI") {
# calculate SE from lower and upper limits
temp_C1 <- (dat$Upr_C[i] - dat$Mean_C[i]) / 1.96
temp_C2 <- (dat$Mean_C[i] - dat$Lwr_C[i]) / 1.96
temp_X1 <- (dat$Upr_X[i] - dat$Mean_X[i]) / 1.96
temp_X2 <- (dat$Mean_X[i] - dat$Lwr_X[i]) / 1.96
# average digitized/recorded values of SE (because we have two) and transform to SD
dat$SD_C[i] <-
mean(abs(c(temp_C1, temp_C2))) * sqrt(dat$N_C[i])
dat$SD_X[i] <-
mean(abs(c(temp_C1, temp_C2))) * sqrt(dat$N_X[i])
} else {
# approximate SD if IQ range is reported
if (dat$Variation.Type[i] == "IQ") {
dat$SD_C[i] <- (dat$Upr_C[i] - dat$Lwr_C[i]) / 1.35
dat$SD_X[i] <- (dat$Upr_X[i] - dat$Lwr_X[i]) / 1.35
}
else {
# if SD is reported leave as such
if (dat$Variation.Type[i] == "SD") {
dat$SD_C[i] <- dat$Variation_C[i]
dat$SD_X[i] <- dat$Variation_X[i]
} else {
# if there is no appropriate data, make these columns NA
dat$SD_C[i] <- NA
dat$SD_X[i] <- NA
}
}
}
}
}Get effect sizes and variances
# within groups standard deviation
dat$S_within <-
sqrt(((dat$N_C - 1) * dat$SD_C ^ 2 + (dat$N_X - 1) * dat$SD_X ^ 2) / (dat$N_C + dat$N_X - 2))
# standardized effect size
for (i in 1:N) {
# standardized mean difference
# we can't include data without variance
if (is.na(dat$S_within[i]) == FALSE & dat$S_within[i] > 0) {
dat$d[i] <- (dat$Mean_X[i] - dat$Mean_C[i]) / dat$S_within[i]
dat$var.d[i] <-
(dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
(2 * (dat$N_C[i] + dat$N_X[i]))
}
else{
# if only F statistic is reported - no cases of this yet!
if (dat$Variation.Type[i] == "F_X") {
dat$d[i] <-
sqrt(dat$F_X[i] * (dat$N_X[i] + dat$N_C[i]) / (dat$N_C[i] * dat$N_X[i]))
dat$var.d[i] <-
(dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
(2 * (dat$N_C[i] + dat$N_X[i]))
} else{
# if Z and N are reported for a regression analysis
if (dat$Variation.Type[i] == "Z_reg") {
r <- dat$Z_X[i] / sqrt(dat$Z_N[i])
dat$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
Vr <- (1 - r ^ 2) ^ 2 / (dat$Z_N[i] - 1)
dat$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
}
else{
# if t and N are reported for a regression analysis
if (dat$Variation.Type[i] == "t") {
r <- sqrt((dat$Z_X[i] ^ 2) / (dat$Z_X[i] ^ 2 + dat$Z_df[i]))
dat$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
Vr <- (1 - r ^ 2) ^ 2 / (dat$R_N[i] - 1)
dat$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
}
else{
# If Z is reported from a comparison of two independent groups
if (dat$Variation.Type[i] == "Z") {
dat$d[i] <-
sqrt(abs(dat$Z_X[i]) * sqrt(dat$N_C[i] + dat$N_X[i]) / (1 - sqrt(
dat$Z_X[i] ^ 2 * (dat$N_C[i] + dat$N_X[i]) ^ -1
)))
dat$var.d[i] <-
(dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
(2 * (dat$N_C[i] + dat$N_X[i]))
}
else{
# converting log OR
if (is.na(dat$Log.OR.v[i]) == FALSE &
dat$Log.OR.v[i] > 0) {
dat$d[i] <- dat$Log.OR[i] * (sqrt(3) / pi)
dat$var.d[i] <- dat$Log.OR.v[i] * 3 / pi ^ 2
} else {
# for now, leave other types of effects as NA
dat$d[i] <- NA
dat$var.d[i] <- NA
}
}
}
}
}
}
}Correct for sample size and get g
# sample size correction factor
for (i in 1:N) {
# if odds ratio or comparison between means
if (is.na(dat$N_C[i]) == F) {
dat$J[i] <- 1 - 3 / (4 * (dat$N_C[i] + dat$N_X[i] - 2) - 1)
} else {
# if correlation
if (is.na(dat$R_N[i]) == F) {
dat$J[i] <- 1 - 3 / (4 * (dat$R_N[i] - 2) - 1)
}else {
dat$J[i] <- NA
}
}
}
# corrected effect size
dat$g <- dat$J * dat$d
# and variance
dat$var.g <- (dat$J ^ 2) * dat$var.dFlip effect sizes for remaining mortality effects
mort <- grep(pattern = "[M,m]ort", x = dat$Trait)
for(i in mort){
dat$g[i] <- -dat$g[i]
}Prepare data
Filter experiments with both epidemiological and infected demographic effects
# For now, get rid of NAs
dat <- dat[(is.na(dat$d) == F),]
# How many experiments do we have?
Exps <- unique(dat$Experiment)
# A data frame of experiments that included Infected demographic and Uninfected demographic effects
dat_1 <- data.frame()
# A data frame of experiments that included Infected demographic and Epidemiological effects
dat_2 <- data.frame()
# populate data sets with corresponding experiments
for (i in Exps) {
tmp <- dat[which(dat$Experiment == i), ]
if ("Uninfected demographic" %in% tmp$Trait.category &
"Infected demographic" %in% tmp$Trait.category) {
dat_1 <- rbind(dat_1, tmp)
}
if ("Epidemiological" %in% tmp$Trait.category &
"Infected demographic" %in% tmp$Trait.category) {
dat_2 <- rbind(dat_2, tmp)
}
}
dat_1 <- dat_1[which(dat_1$Trait.category == "Uninfected demographic" | dat_1$Trait.category == "Infected demographic"),]
dat_2 <- dat_2[which(dat_2$Trait.category == "Epidemiological" | dat_2$Trait.category == "Infected demographic"),]
dat_2$Trait.type <- factor(dat_2$Trait.type, levels = c("Prevalence", "Intensity", "Fecundity", "Survivorship"), ordered = T)Variance-covariance of sampling errors matrix (for multiple comparisons)
This code produces a variance-covariance of sampling errors matrix, n x n, with n = number of effect sizes. It requires the dataset (dat), Hedge’s g (g), and the variance of Hedge’s g (var.g) to be defined above. It uses “Experiment”, “Level_C”, “Trait”, “N_C”, and “N_X” columns.
Identify negative eigenvalues
First, we find out which studies have negative eigenvalues in their variance-covar matrix. This means covariance between treatment effects with the same control is larger than variance and it’s an indicator of suspiciously low variance and potentially an incorrectly reported N or an error while digitizing data.
For data set 1
studies_1 <- unique(dat_1$Study)
studies_to_check <- c()
negative_eigen <- c()
for (k in studies_1) {
dat_study <- dat_1[which(dat_1$Study == k),]
varcovmat = matrix(0,
nrow = dim(dat_study)[1],
ncol = dim(dat_study)[1])
for (i in 1:dim(dat_study)[1]) {
for (j in 1:dim(dat_study)[1]) {
if (i == j) {
varcovmat[i, j] = dat_study$var.g[i]
} else{
if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
dat_study[i, "Level_C"] == dat_study[j, "Level_C"] &
dat_study[i, "Trait"] == dat_study[j, "Trait"] &
dat_study[i, "Trait.category"] == dat_study[j, "Trait.category"]) {
varcovmat[i, j] = 1 / dat_study[i, "N_C"] + dat_study$g[i] * dat_study$g[j] /
(dat_study[i, "N_C"] + dat_study[i, "N_X"] + dat_study[j, "N_X"])
}
}
}
}
val <- eigen(varcovmat)
for (m in 1:length(val$values)) {
if (val$values[m] < 0) {
studies_to_check <- append(studies_to_check, k)
negative_eigen <- append(negative_eigen, val$values[m])
}
}
}
Neg_eigen_1 <- data.frame(study = studies_to_check, eigen = negative_eigen)
Neg_eigen_1 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| study | eigen |
|---|---|
| Shostak et al. 2015. Journal of Parasitology | -1.0728652 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0311074 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0459004 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0506383 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.1843072 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.2369469 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3778444 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.5191928 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.5801272 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.6831878 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7934325 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.8845345 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9065495 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9289586 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9316950 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9469193 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.0721254 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.1539142 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.2820301 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.3782808 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.4589050 |
| Schluter-Vorberg et al. 2017. Aquatic Toxicology | -0.0496863 |
| Schluter-Vorberg et al. 2017. Aquatic Toxicology | -0.0753315 |
| Waki & Yoshinaga 2015. FishPathology | -0.1246226 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.0097697 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.0487760 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.1770843 |
| Pochini and Hoverman. 2017. Environmental Pollution | -0.1229350 |
| Pochini and Hoverman. 2017. Environmental Pollution | -1.1461176 |
| Civitello et al. 2020. Proc B | -0.5611703 |
| Civitello et al. 2020. Proc B | -3.8834534 |
For data set 2
studies_2 <- unique(dat_2$Study)
studies_to_check <- c()
negative_eigen <- c()
for (k in studies_2) {
dat_study <- dat_2[which(dat_2$Study == k),]
varcovmat = matrix(0,
nrow = dim(dat_study)[1],
ncol = dim(dat_study)[1])
for (i in 1:dim(dat_study)[1]) {
for (j in 1:dim(dat_study)[1]) {
if (i == j) {
varcovmat[i, j] = dat_study$var.g[i]
} else{
if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
dat_study[i, "Level_C"] == dat_study[j, "Level_C"] &
dat_study[i, "Trait"] == dat_study[j, "Trait"] &
dat_study[i, "Trait.category"] == dat_study[j, "Trait.category"]) {
varcovmat[i, j] = 1 / dat_study[i, "N_C"] + dat_study$g[i] * dat_study$g[j] /
(dat_study[i, "N_C"] + dat_study[i, "N_X"] + dat_study[j, "N_X"])
}
}
}
}
val <- eigen(varcovmat)
for (m in 1:length(val$values)) {
if (val$values[m] < 0) {
studies_to_check <- append(studies_to_check, k)
negative_eigen <- append(negative_eigen, val$values[m])
}
}
}
Neg_eigen_2 <- data.frame(study = studies_to_check, eigen = negative_eigen)
Neg_eigen_2 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| study | eigen |
|---|---|
| Bradley et al. 2019. Plos One | -0.0015336 |
| Bradley et al. 2019. Plos One | -0.0250523 |
| Bradley et al. 2019. Plos One | -0.0379281 |
| Bradley et al. 2019. Plos One | -0.1081902 |
| Souto et al. 2015. VetMicrobiol. | -0.5585385 |
| Carrington et al. 2013. PLOS NegTropDiseases | -0.0050917 |
| Carrington et al. 2013. PLOS NegTropDiseases | -0.0603175 |
| Carrington et al. 2013. PLOS NegTropDiseases | -0.1012185 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0311074 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0506383 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.2369469 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.2917917 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3005591 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3576572 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3844600 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3890634 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3998072 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.5191928 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.5801272 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.6662958 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.6831878 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7429232 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7522409 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7660915 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7758519 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.7934325 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.8008515 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9216851 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9233991 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9253113 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9316950 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.9398650 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.0283485 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.0694983 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.0721254 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.1094691 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.2178041 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.3472667 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.3815538 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.4370017 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.5357111 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.5532399 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.7623978 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -1.9073762 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -2.3818336 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -2.6483026 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -6.8202809 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -7.0495394 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -8.3833406 |
| Schluter-Vorberg et al. 2017. Aquatic Toxicology | -0.0496863 |
| Schluter-Vorberg et al. 2017. Aquatic Toxicology | -0.0753315 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.0077688 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.0487760 |
| Waki&Yoshinaga. 2013. Fish Sci | -0.1770843 |
| Purcell et al. 2016 J Fish Diseases | -0.1184091 |
Based on these eigenvalues, we looked back at the studies/experiments and determined if they had issues with reporting N or variance. We decided to remove two such studies for data set 1 (Ashraf et al. 2017 and Shostak et al. 2015) and one such study for data set 2 (Ashraf et al. 2017).
As we found no errors or causes of concern for other studies with barely negative eigenvalues, we are forcing the var-covar matrices to the nearest positive definite values, while keeping the diagonals (var) to the estimated values and modifying only the expected covariances.
Also for now we are excluding Experiment 121 from Civitello et al 2020 Proc. B. This study has a relatively low variance for the control of the uninfected demographic fecundity effects, which results in a negative eigenvalue.
# exclude suspicious studies with highly negative eigenvalues
dat_1 <- dat_1[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_1$Study),]
dat_1 <- dat_1[-grep("Shostak et al. 2015. Journal of Parasitology", dat_1$Study),]
dat_1 <- dat_1[-grep(121, dat_1$Experiment),]
varcovmat_1 = matrix(0, nrow = dim(dat_1)[1], ncol = dim(dat_1)[1])
for (i in 1:dim(dat_1)[1]) {
for (j in 1:dim(dat_1)[1]) {
if (i == j) {varcovmat_1[i,j] = dat_1$var.g[i]}else{
if (dat_1[i, "Experiment"] == dat_1[j, "Experiment"] & dat_1[i, "Level_C"] == dat_1[j, "Level_C"] & dat_1[i, "Trait"] == dat_1[j, "Trait"] & dat_1[i, "Trait.category"] == dat_1[j, "Trait.category"]) {
varcovmat_1[i,j] = 1/dat_1[i,"N_C"] + dat_1$g[i]*dat_1$g[j]/(dat_1[i,"N_C"] + dat_1[i,"N_X"] + dat_1[j,"N_X"])
}
}
}
}
# correct negative eigens in the few studies with large covar to var ratios
varcovmat_1_PD <- nearPD(varcovmat_1, keepDiag = TRUE)
# exclude suspicious studies with highly negative eigenvalues
dat_2 <- dat_2[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_2$Study),]
# Carrington et al. 2013 has extremely large sampling variance for a parasite prevalence effect. DOUBLE CHECK!!!
dat_2 <- dat_2[-grep("Carrington et al. 2013. PLOS NegTropDiseases", dat_2$Study),]
varcovmat_2 = matrix(0, nrow = dim(dat_2)[1], ncol = dim(dat_2)[1])
for (i in 1:dim(dat_2)[1]) {
for (j in 1:dim(dat_2)[1]) {
if (i == j) {varcovmat_2[i,j] = dat_2$var.g[i]}else{
if (dat_2[i, "Experiment"] == dat_2[j, "Experiment"] & dat_2[i, "Level_C"] == dat_2[j, "Level_C"] & dat_2[i, "Trait"] == dat_2[j, "Trait"] & dat_2[i, "Trait.category"] == dat_2[j, "Trait.category"]) {
varcovmat_2[i,j] = 1/dat_2[i,"N_C"] + dat_2$g[i]*dat_2$g[j]/(dat_2[i,"N_C"] + dat_2[i,"N_X"] + dat_2[j,"N_X"])
}
}
}
}
# correct eigenvalues in the few studies with large covar to var ratios
varcovmat_2_PD <- nearPD(varcovmat_2, keepDiag = TRUE)Create a new column with a category of vertebrate or invertebrate (keeping this for now)
dat_1 <- mutate(dat_1, Host.type.2 = case_when(
Host.type == "Fish" ~ "Vertebrate",
Host.type == "Arthropod" ~ "Invertebrate",
Host.type == "Amphibian" ~ "Vertebrate",
Host.type == "Mollusc" ~ "Invertebrate"))
dat_2 <- mutate(dat_2, Host.type.2 = case_when(
Host.type == "Fish" ~ "Vertebrate",
Host.type == "Arthropod" ~ "Invertebrate",
Host.type == "Amphibian" ~ "Vertebrate",
Host.type == "Mollusc" ~ "Invertebrate",
Host.type == "Reptile" ~ "Vertebrate",
Host.type == "Bird" ~ "Vertebrate",
Host.type == "Mammal" ~ "Vertebrate"))Numbers mentioned in the methods
The Methods section describes various features of the data, indicating the number of studies/effects with those features. Numbers are calculated here.
# First, get all the curated data
final_dat <- merge(dat_2, dat_1, all = T )
nrow(final_dat)## [1] 891# do we have any observational studies left after filtering?
levels(as.factor(final_dat$Study.type))## [1] "Experimental"# Number of studies that report median instead of mean
c(unique(non_par$Study), unique(non_par_2$Study)) %in% final_dat$Study## [1] TRUE TRUE TRUE TRUElength(unique(non_par$Study)) + length(unique(non_par_2$Study)) ## [1] 4length(non_par$ID) + length(non_par_2$ID) ## [1] 13# number of studies with only data range
length(non_par$ID)## [1] 8length(unique(non_par$Study))## [1] 1# number of studies with only IQ
length(non_par_2$ID)## [1] 5length(unique(non_par_2$Study))## [1] 3# number of studies with OR
length(unique(final_dat$Study[final_dat$Variation.Type == "OR"]))## [1] 66# number or studies w/ multiple treatment levels or measuring time points
k = c()
for(i in 2:nrow(final_dat)){
if(final_dat$Trait[i] == final_dat$Trait[i-1] &
final_dat$X_gradient[i] == final_dat$X_gradient[i-1] &
final_dat$Host[i] == final_dat$Host[i-1] &
final_dat$Parasite[i] == final_dat$Parasite[i-1]) {
k <- append(k, final_dat$Experiment[i])
}
}
length(unique(k))## [1] 110# studies with negative eigen values
sum(unique(c(Neg_eigen_1$eigen,Neg_eigen_2$study)) %in% final_dat$Study == T)## [1] 5# how many experiments
length(unique(final_dat$Experiment))## [1] 122# how many studies include more than one experiment
k = c()
for(i in 2:nrow(final_dat)){
if(final_dat$Experiment[i] != final_dat$Experiment[i-1] &
final_dat$Study[i] == final_dat$Study[i-1]) {
k <- append(k, final_dat$Study[i])
}
}
length(unique(k))## [1] 21# how many host species per taxonomic groups
Host_count <- final_dat %>% group_by(Host.type) %>%
summarise(N_taxa = length(unique(Host)))
Host_count %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Host.type | N_taxa |
|---|---|
| Amphibian | 21 |
| Arthropod | 20 |
| Bird | 2 |
| Fish | 13 |
| Mammal | 1 |
| Mollusc | 13 |
| Reptile | 1 |
# how many parasite taxa per taxonomic groups
Par_count <- final_dat %>% group_by(Parasite.type) %>%
summarise(N_taxa = length(unique(Parasite)))
Par_count %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Parasite.type | N_taxa |
|---|---|
| Bacteria | 14 |
| Fungus | 6 |
| Helminth | 12 |
| Myxozoa | 1 |
| Protozoan | 8 |
| Virus | 23 |
Analysis
Question 1 - Is there a differencial response to stressors between infected and uninfected host?
Model selection
# First, we evaluate some code that generates helper functions needed so that metafor and MuMIn could interact as necessary
eval(metafor:::.MuMIn)
# Now, we take the full model and fit the rest of the models and examine those models whose AICc value is no more than 10 units away from that of the best model
full_model.1 <- rma.mv(
g,
V = varcovmat_1_PD$mat,
mods = ~ Trait.category*Trait.type*Gradient.category,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "ML"
)
model_selection.1 <- dredge(full_model.1, trace = 2) ## Fixed term is "(Intercept)"##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|====== | 9%
|
|======== | 12%
|
|=========== | 16%
|
|============= | 18%
|
|================= | 24%
|
|===================== | 30%
|
|========================== | 37%
|
|============================== | 43%
|
|================================== | 49%
|
|===================================================================== | 99%subset(model_selection.1, delta <= 10, recalc.weights = FALSE)## Global model call: rma.mv(yi = g, V = varcovmat_1_PD$mat, mods = ~Trait.category *
## Trait.type * Gradient.category, random = list(~1 | Experiment,
## ~1 | ID), data = dat_1, method = "ML")
## ---
## Model selection table
## (Int) Grd.ctg Trt.ctg Trt.typ Grd.ctg:Trt.ctg Grd.ctg:Trt.typ
## 22 + + + +
## 24 + + + + +
## 56 + + + + +
## 1 +
## 32 + + + + + +
## 5 + +
## 3 + +
## 7 + + +
## 64 + + + + + +
## 2 + +
## 6 + + +
## 4 + + +
## 39 + + +
## 8 + + + +
## 128 + + + + + +
## Trt.ctg:Trt.typ Grd.ctg:Trt.ctg:Trt.typ df logLik AICc delta weight
## 22 8 -555.340 1127.1 0.00 0.405
## 24 9 -554.802 1128.1 1.02 0.243
## 56 + 10 -554.791 1130.2 3.11 0.086
## 1 3 -562.492 1131.0 3.98 0.055
## 32 11 -554.509 1131.7 4.66 0.039
## 5 4 -561.842 1131.8 4.72 0.038
## 3 4 -561.889 1131.9 4.82 0.036
## 7 5 -561.316 1132.8 5.73 0.023
## 64 + 12 -554.457 1133.8 6.69 0.014
## 2 5 -561.802 1133.8 6.70 0.014
## 6 6 -561.026 1134.3 7.21 0.011
## 4 6 -561.168 1134.6 7.49 0.010
## 39 + 6 -561.313 1134.8 7.78 0.008
## 8 7 -560.480 1135.3 8.19 0.007
## 128 + + 14 -553.535 1136.2 9.14 0.004
## Models ranked by AICc(x)# relative importance values for the predictors can be obtained with:
sw(model_selection.1)## Trait.type Gradient.category Gradient.category:Trait.type
## Sum of weights: 0.88 0.84 0.79
## N containing models: 14 14 6
## Trait.category Trait.category:Trait.type
## Sum of weights: 0.48 0.12
## N containing models: 14 6
## Gradient.category:Trait.category
## Sum of weights: 0.06
## N containing models: 6
## Gradient.category:Trait.category:Trait.type
## Sum of weights: <0.01
## N containing models: 1Full Model
We fit the full model to show there is no effect of infection status
Q1.f <-
rma.mv(
g ~ Trait.category:Trait.type:Gradient.category -1,
V = varcovmat_1_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
pal <- RColorBrewer::brewer.pal(9,"Spectral")
Q1.f_OP <-
orchard_plot(
Q1.f,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
scale_fill_manual(values = rep(pal, each = 2)) + scale_colour_manual(values = rep(pal, each = 2)) +
theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_y_discrete(
labels = c(
"Infected:EE:fecund",
"Uninfected:EE:fecund",
"Infected:EE:surv",
"Uninfected:EE:surv",
"Infected:CP:fecund",
"Uninfected:CP:fecund",
"Infected:CP:surv",
"Uninfected:CP:surv",
"Infected:RL:fecund",
"Uninfected:RL:fecund",
"Infected:RL:surv",
"Uninfected:RL:surv"
)
)
Q1.f_OP
Best Model: Effects in survival vs fecundity under different types of stress, combining infected and uninfected
Q1 <-
rma.mv(
g ~ Trait.type:Gradient.category -1 ,
V = varcovmat_1_PD$mat,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
summary(Q1)##
## Multivariate Meta-Analysis Model (k = 384; method: REML)
##
## logLik Deviance AIC BIC AICc
## -548.7330 1097.4659 1113.4659 1144.9451 1113.8562
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2560 0.5060 78 no Experiment
## sigma^2.2 0.3607 0.6005 384 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 378) = 217601.2265, p-val < .0001
##
## Test of Moderators (coefficients 1:6):
## QM(df = 6) = 46.7815, p-val < .0001
##
## Model Results:
##
## estimate se zval
## Trait.typeFecundity:Gradient.categoryEnvironment -0.3255 0.2429 -1.3404
## Trait.typeSurvivorship:Gradient.categoryEnvironment -0.5161 0.1016 -5.0776
## Trait.typeFecundity:Gradient.categoryPollution -0.1406 0.3153 -0.4458
## Trait.typeSurvivorship:Gradient.categoryPollution -0.3510 0.1477 -2.3767
## Trait.typeFecundity:Gradient.categoryResource -0.8637 0.2161 -3.9974
## Trait.typeSurvivorship:Gradient.categoryResource 0.0025 0.2050 0.0121
## pval ci.lb ci.ub
## Trait.typeFecundity:Gradient.categoryEnvironment 0.1801 -0.8015 0.1505
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 -0.7153 -0.3169
## Trait.typeFecundity:Gradient.categoryPollution 0.6557 -0.7586 0.4775
## Trait.typeSurvivorship:Gradient.categoryPollution 0.0175 -0.6404 -0.0615
## Trait.typeFecundity:Gradient.categoryResource <.0001 -1.2872 -0.4402
## Trait.typeSurvivorship:Gradient.categoryResource 0.9904 -0.3994 0.4044
##
## Trait.typeFecundity:Gradient.categoryEnvironment
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryPollution
## Trait.typeSurvivorship:Gradient.categoryPollution *
## Trait.typeFecundity:Gradient.categoryResource ***
## Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q1$beta, sei = Q1$se, ci.lb = Q1$ci.lb, ci.ub = Q1$ci.ub,
annotate=TRUE, showweights=T, header=F,
slab = c("Environment fecundity", "Environment survivorship", "Pollution fecundity", "Pollution survivorship", "Resource fecundity", "Resource survivorship"))
Make Table 2 with results from Q1
Table2 <- data.frame(Stressor_type = c("Endogenous environment", "Endogenous environment", "Chemical pollution", "Chemical pollution", "Resource limitation","Resource limitation"),
Response_trait = c(rep(c("Fecundity", "Survivorship"), 3)),
Overall_mean = round(Q1$beta,3),
Lower_95 = round(Q1$ci.lb,3),
Upper_95 = round(Q1$ci.ub,3),
P_value = round(Q1$pval,3),
N_effects = c(length(dat_1$ID[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"]),
length(dat_1$ID[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"]),
length(dat_1$ID[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"]),
length(dat_1$ID[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"]),
length(dat_1$ID[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"]),
length(dat_1$ID[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"])),
N_experiments = c(length(unique(dat_1$Experiment[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Experiment[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"])),
length(unique(dat_1$Experiment[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Experiment[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"])),
length(unique(dat_1$Experiment[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Experiment[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"]))
),
N_host_taxa = c(length(unique(dat_1$Host[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Host[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"])),
length(unique(dat_1$Host[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Host[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"])),
length(unique(dat_1$Host[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"])),
length(unique(dat_1$Host[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"]))
))
rownames(Table2) <- NULL
Table2 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | N_effects | N_experiments | N_host_taxa |
|---|---|---|---|---|---|---|---|---|
| Endogenous environment | Fecundity | -0.326 | -0.802 | 0.150 | 0.180 | 14 | 3 | 3 |
| Endogenous environment | Survivorship | -0.516 | -0.715 | -0.317 | 0.000 | 189 | 45 | 34 |
| Chemical pollution | Fecundity | -0.141 | -0.759 | 0.477 | 0.656 | 14 | 4 | 2 |
| Chemical pollution | Survivorship | -0.351 | -0.640 | -0.062 | 0.017 | 76 | 23 | 15 |
| Resource limitation | Fecundity | -0.864 | -1.287 | -0.440 | 0.000 | 59 | 6 | 4 |
| Resource limitation | Survivorship | 0.002 | -0.399 | 0.404 | 0.990 | 32 | 6 | 5 |
Make Figure 2 with the results of the best model for Q1
pal <- RColorBrewer::brewer.pal(11,"Spectral")
Q1_OP <-
orchard_plot(
Q1,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
scale_fill_manual(values = rep(pal[c(8,11,5)], each = 2)) + scale_colour_manual(values = rep(pal[c(8,11,5)], each = 2)) +
theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_y_discrete(
labels = c("EE:fecundity",
"EE:survivorship",
"CP:fecundity",
"CP:survivorship",
"RL:fecundity",
"RL:survivorship"))## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.Q1_OP
RVE
We check if our results are qualitatively similar if we use the robust variance estimator instead of the var-covar matrix for effects estimated in the same experiment.
Q1.R <-
rma.mv(
g ~ Trait.type:Gradient.category -1 ,
V = dat_1$var.g,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
Q1.R_est <- clubSandwich::coef_test(Q1.R, vcov = "CR2")## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwichforest.default(x= Q1.R_est$beta, sei = Q1.R_est$SE, pch = 19,
annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"),
slab = c("EE:fecundity", "EE:survivorship", "CP:fecundity", "CP:survivorship", "RL:fecundity", "RL:survivorship"),
xlab = "Standardised mean difference (Hedge's g)")
Question 2 - Do fitness and infectivity traits have different sensitivity to stressors?
Model selection
All stressors
full_model.2 <- rma.mv(
g,
V = varcovmat_2_PD$mat,
mods = ~ Trait.type*Gradient.category,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "ML"
)
# Define if using Trait.type or Trait.category
model_selection.2 <- dredge(full_model.2, trace = 2) ## Fixed term is "(Intercept)"##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|============================================================= | 88%subset(model_selection.2, delta <= 10, recalc.weights = FALSE)## Global model call: rma.mv(yi = g, V = varcovmat_2_PD$mat, mods = ~Trait.type * Gradient.category,
## random = list(~1 | Experiment, ~1 | ID), data = dat_2, method = "ML")
## ---
## Model selection table
## (Int) Grd.ctg Trt.typ Grd.ctg:Trt.typ df logLik AICc delta weight
## 8 + + + + 14 -1042.743 2114.1 0 1
## Models ranked by AICc(x)# relative importance values for the predictors can be obtained with:
sw(model_selection.2)## Trait.type Gradient.category Gradient.category:Trait.type
## Sum of weights: 1 1 1
## N containing models: 3 3 1Fitness versus infectivity traits
full_model.2.rvf <- rma.mv(
g,
V = varcovmat_2_PD$mat,
mods = ~ Trait.category*Gradient.category,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "ML"
)
# Define if using Trait.type or Trait.category
model_selection.2.rvf <- dredge(full_model.2.rvf, trace = 2) ## Fixed term is "(Intercept)"##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|============================================================= | 88%subset(model_selection.2.rvf, delta <= 10, recalc.weights = FALSE)## Global model call: rma.mv(yi = g, V = varcovmat_2_PD$mat, mods = ~Trait.category *
## Gradient.category, random = list(~1 | Experiment, ~1 | ID),
## data = dat_2, method = "ML")
## ---
## Model selection table
## (Int) Grd.ctg Trt.ctg Grd.ctg:Trt.ctg df logLik AICc delta weight
## 8 + + + + 8 -1053.678 2123.6 0 1
## Models ranked by AICc(x)# relative importance values for the predictors can be obtained with:
sw(model_selection.2.rvf)## Trait.category Gradient.category
## Sum of weights: 1 1
## N containing models: 3 3
## Gradient.category:Trait.category
## Sum of weights: 1
## N containing models: 1Best model: Trait type x Gradient category interaction (lowest AICc)
Q2 <-
rma.mv(
g ~ Trait.type:Gradient.category -1,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
summary(Q2)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -1031.4252 2062.8504 2090.8504 2154.0356 2091.4877
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1697 0.4120 113 no Experiment
## sigma^2.2 0.4486 0.6698 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 674) = 1423699.2554, p-val < .0001
##
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 103.7416, p-val < .0001
##
## Model Results:
##
## estimate se zval
## Trait.typePrevalence:Gradient.categoryEnvironment 0.1010 0.1293 0.7813
## Trait.typeIntensity:Gradient.categoryEnvironment 0.3978 0.0937 4.2434
## Trait.typeFecundity:Gradient.categoryEnvironment 0.0318 0.2802 0.1134
## Trait.typeSurvivorship:Gradient.categoryEnvironment -0.5289 0.0901 -5.8734
## Trait.typePrevalence:Gradient.categoryPollution -0.5227 0.2520 -2.0741
## Trait.typeIntensity:Gradient.categoryPollution -0.0192 0.1266 -0.1518
## Trait.typeFecundity:Gradient.categoryPollution 0.0887 0.2668 0.3324
## Trait.typeSurvivorship:Gradient.categoryPollution -0.2627 0.1318 -1.9933
## Trait.typePrevalence:Gradient.categoryResource -0.1503 0.2465 -0.6096
## Trait.typeIntensity:Gradient.categoryResource -0.4346 0.1826 -2.3803
## Trait.typeFecundity:Gradient.categoryResource -0.7416 0.2021 -3.6699
## Trait.typeSurvivorship:Gradient.categoryResource 0.0049 0.1957 0.0253
## pval ci.lb ci.ub
## Trait.typePrevalence:Gradient.categoryEnvironment 0.4347 -0.1524 0.3545
## Trait.typeIntensity:Gradient.categoryEnvironment <.0001 0.2141 0.5815
## Trait.typeFecundity:Gradient.categoryEnvironment 0.9097 -0.5173 0.5809
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 -0.7054 -0.3524
## Trait.typePrevalence:Gradient.categoryPollution 0.0381 -1.0166 -0.0288
## Trait.typeIntensity:Gradient.categoryPollution 0.8793 -0.2673 0.2289
## Trait.typeFecundity:Gradient.categoryPollution 0.7396 -0.4343 0.6116
## Trait.typeSurvivorship:Gradient.categoryPollution 0.0462 -0.5210 -0.0044
## Trait.typePrevalence:Gradient.categoryResource 0.5421 -0.6335 0.3329
## Trait.typeIntensity:Gradient.categoryResource 0.0173 -0.7924 -0.0767
## Trait.typeFecundity:Gradient.categoryResource 0.0002 -1.1376 -0.3455
## Trait.typeSurvivorship:Gradient.categoryResource 0.9799 -0.3786 0.3885
##
## Trait.typePrevalence:Gradient.categoryEnvironment
## Trait.typeIntensity:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryEnvironment
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typePrevalence:Gradient.categoryPollution *
## Trait.typeIntensity:Gradient.categoryPollution
## Trait.typeFecundity:Gradient.categoryPollution
## Trait.typeSurvivorship:Gradient.categoryPollution *
## Trait.typePrevalence:Gradient.categoryResource
## Trait.typeIntensity:Gradient.categoryResource *
## Trait.typeFecundity:Gradient.categoryResource ***
## Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q2$beta, sei = Q2$se, ci.lb = Q2$ci.lb, ci.ub = Q2$ci.ub,
annotate=TRUE, showweights=T, header=F,
slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
"CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
"RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))
par_short <- c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
"CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
"RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship")
pal <- RColorBrewer::brewer.pal(11,"Spectral")
Q2_OP <-
orchard_plot(
Q2,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
theme(legend.position = c(0.3, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_fill_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 2)) +
scale_colour_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 2)) +
scale_y_discrete(labels = par_short)## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.Q2_OP
Contrasts
Flipped values for fitness and infectivity Hedge’s g
dat_2$flipped <- dat_2$g
flip <- grep(pattern = "Intensity|Prevalence", x = dat_2$Trait.type)
for(i in flip){
dat_2$flipped[i] <- -dat_2$flipped[i]
}Re-run model with flipped values
Q2.flip <-
rma.mv(
flipped ~ Trait.type:Gradient.category -1,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
summary(Q2.flip)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -979.7187 1959.4374 1987.4374 2050.6227 1988.0748
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4654 0.6822 113 no Experiment
## sigma^2.2 0.3018 0.5493 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 674) = 1423699.2554, p-val < .0001
##
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 67.2068, p-val < .0001
##
## Model Results:
##
## estimate se zval
## Trait.typePrevalence:Gradient.categoryEnvironment -0.3614 0.1363 -2.6525
## Trait.typeIntensity:Gradient.categoryEnvironment -0.4498 0.1086 -4.1415
## Trait.typeFecundity:Gradient.categoryEnvironment -0.4919 0.2574 -1.9108
## Trait.typeSurvivorship:Gradient.categoryEnvironment -0.6301 0.1047 -6.0168
## Trait.typePrevalence:Gradient.categoryPollution 0.3523 0.2609 1.3503
## Trait.typeIntensity:Gradient.categoryPollution 0.0600 0.1502 0.3995
## Trait.typeFecundity:Gradient.categoryPollution -0.2026 0.2611 -0.7759
## Trait.typeSurvivorship:Gradient.categoryPollution -0.2314 0.1525 -1.5172
## Trait.typePrevalence:Gradient.categoryResource 0.1854 0.2365 0.7840
## Trait.typeIntensity:Gradient.categoryResource 0.2959 0.1908 1.5506
## Trait.typeFecundity:Gradient.categoryResource -0.5041 0.2108 -2.3912
## Trait.typeSurvivorship:Gradient.categoryResource -0.2563 0.1940 -1.3214
## pval ci.lb ci.ub
## Trait.typePrevalence:Gradient.categoryEnvironment 0.0080 -0.6285 -0.0944
## Trait.typeIntensity:Gradient.categoryEnvironment <.0001 -0.6626 -0.2369
## Trait.typeFecundity:Gradient.categoryEnvironment 0.0560 -0.9964 0.0127
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 -0.8354 -0.4248
## Trait.typePrevalence:Gradient.categoryPollution 0.1769 -0.1591 0.8637
## Trait.typeIntensity:Gradient.categoryPollution 0.6896 -0.2344 0.3543
## Trait.typeFecundity:Gradient.categoryPollution 0.4378 -0.7145 0.3092
## Trait.typeSurvivorship:Gradient.categoryPollution 0.1292 -0.5303 0.0675
## Trait.typePrevalence:Gradient.categoryResource 0.4331 -0.2782 0.6490
## Trait.typeIntensity:Gradient.categoryResource 0.1210 -0.0781 0.6699
## Trait.typeFecundity:Gradient.categoryResource 0.0168 -0.9173 -0.0909
## Trait.typeSurvivorship:Gradient.categoryResource 0.1864 -0.6365 0.1239
##
## Trait.typePrevalence:Gradient.categoryEnvironment **
## Trait.typeIntensity:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryEnvironment .
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typePrevalence:Gradient.categoryPollution
## Trait.typeIntensity:Gradient.categoryPollution
## Trait.typeFecundity:Gradient.categoryPollution
## Trait.typeSurvivorship:Gradient.categoryPollution
## Trait.typePrevalence:Gradient.categoryResource
## Trait.typeIntensity:Gradient.categoryResource
## Trait.typeFecundity:Gradient.categoryResource *
## Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q2.flip$beta, sei = Q2.flip$se, ci.lb = Q2.flip$ci.lb, ci.ub = Q2.flip$ci.ub,
annotate=TRUE, showweights=T, header=F,
slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
"CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
"RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))
Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts with flipped values.
# Environment
# Prevalence vs Fecundity
anova(Q2.flip, X=c(-1,0,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.1304 0.2672 -0.4882 0.6254
##
## Test of Hypothesis:
## QM(df = 1) = 0.2383, p-val = 0.6254# Prevalence vs Survivorship
anova(Q2.flip, X=c(-1,0,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.2687 0.1231 -2.1823 0.0291
##
## Test of Hypothesis:
## QM(df = 1) = 4.7626, p-val = 0.0291# Intensity vs Fecundity
anova(Q2.flip, X=c(0,-1,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.0421 0.2567 -0.1641 0.8697
##
## Test of Hypothesis:
## QM(df = 1) = 0.0269, p-val = 0.8697# Intensity vs Survivorship
anova(Q2.flip, X=c(0,-1,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.1803 0.0986 -1.8294 0.0673
##
## Test of Hypothesis:
## QM(df = 1) = 3.3468, p-val = 0.0673# Pollution
# Prevalence vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,-1,0,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.5549 0.3376 -1.6437 0.1002
##
## Test of Hypothesis:
## QM(df = 1) = 2.7019, p-val = 0.1002# Prevalence vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,-1,0,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.5837 0.2564 -2.2764 0.0228
##
## Test of Hypothesis:
## QM(df = 1) = 5.1819, p-val = 0.0228# Intensity vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,-1,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.2626 0.2459 -1.0678 0.2856
##
## Test of Hypothesis:
## QM(df = 1) = 1.1402, p-val = 0.2856# Intensity vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,-1,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.2914 0.1328 -2.1938 0.0283
##
## Test of Hypothesis:
## QM(df = 1) = 4.8126, p-val = 0.0283# Resources
# Prevalence vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,-1,0,1,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: -0.6895 0.2555 -2.6991 0.0070
##
## Test of Hypothesis:
## QM(df = 1) = 7.2852, p-val = 0.0070# Prevalence vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,-1,0,0,1))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: -0.4417 0.2596 -1.7013 0.0889
##
## Test of Hypothesis:
## QM(df = 1) = 2.8943, p-val = 0.0889# Intensity vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,0,-1,1,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: -0.8000 0.2025 -3.9509 <.0001
##
## Test of Hypothesis:
## QM(df = 1) = 15.6098, p-val < .0001# Intensity vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,0,-1,0,1))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: -0.5522 0.2128 -2.5950 0.0095
##
## Test of Hypothesis:
## QM(df = 1) = 6.7338, p-val = 0.0095Absolute values for fitness and infectivity Hedge’s g
#absolute <- grep(pattern = "Fecundity|Survivorship|Intensity|Prevalence", x = dat_2$Trait.type)
#for(i in absolute){
# dat_2$g[i] <- abs(dat_2$g[i])
#}
dat_2$absolute <- abs(dat_2$g)Re-run model with absolute values
Q2.abs <-
rma.mv(
absolute ~ Trait.type:Gradient.category -1,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
summary(Q2.abs)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -781.8300 1563.6600 1591.6600 1654.8452 1592.2973
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1852 0.4304 113 no Experiment
## sigma^2.2 0.1073 0.3275 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 674) = 1427007.9678, p-val < .0001
##
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 269.2974, p-val < .0001
##
## Model Results:
##
## estimate se zval
## Trait.typePrevalence:Gradient.categoryEnvironment 0.8172 0.0997 8.1933
## Trait.typeIntensity:Gradient.categoryEnvironment 0.9734 0.0771 12.6211
## Trait.typeFecundity:Gradient.categoryEnvironment 0.5422 0.1850 2.9299
## Trait.typeSurvivorship:Gradient.categoryEnvironment 0.7927 0.0738 10.7403
## Trait.typePrevalence:Gradient.categoryPollution 0.7153 0.2012 3.5558
## Trait.typeIntensity:Gradient.categoryPollution 0.6677 0.1039 6.4237
## Trait.typeFecundity:Gradient.categoryPollution 0.8958 0.1966 4.5568
## Trait.typeSurvivorship:Gradient.categoryPollution 0.6026 0.1060 5.6841
## Trait.typePrevalence:Gradient.categoryResource 0.5347 0.1773 3.0152
## Trait.typeIntensity:Gradient.categoryResource 0.8019 0.1376 5.8282
## Trait.typeFecundity:Gradient.categoryResource 0.8126 0.1521 5.3440
## Trait.typeSurvivorship:Gradient.categoryResource 0.8107 0.1406 5.7651
## pval ci.lb ci.ub
## Trait.typePrevalence:Gradient.categoryEnvironment <.0001 0.6217 1.0127
## Trait.typeIntensity:Gradient.categoryEnvironment <.0001 0.8222 1.1245
## Trait.typeFecundity:Gradient.categoryEnvironment 0.0034 0.1795 0.9048
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 0.6480 0.9374
## Trait.typePrevalence:Gradient.categoryPollution 0.0004 0.3210 1.1096
## Trait.typeIntensity:Gradient.categoryPollution <.0001 0.4640 0.8714
## Trait.typeFecundity:Gradient.categoryPollution <.0001 0.5105 1.2811
## Trait.typeSurvivorship:Gradient.categoryPollution <.0001 0.3948 0.8103
## Trait.typePrevalence:Gradient.categoryResource 0.0026 0.1871 0.8822
## Trait.typeIntensity:Gradient.categoryResource <.0001 0.5322 1.0715
## Trait.typeFecundity:Gradient.categoryResource <.0001 0.5145 1.1106
## Trait.typeSurvivorship:Gradient.categoryResource <.0001 0.5351 1.0863
##
## Trait.typePrevalence:Gradient.categoryEnvironment ***
## Trait.typeIntensity:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryEnvironment **
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typePrevalence:Gradient.categoryPollution ***
## Trait.typeIntensity:Gradient.categoryPollution ***
## Trait.typeFecundity:Gradient.categoryPollution ***
## Trait.typeSurvivorship:Gradient.categoryPollution ***
## Trait.typePrevalence:Gradient.categoryResource **
## Trait.typeIntensity:Gradient.categoryResource ***
## Trait.typeFecundity:Gradient.categoryResource ***
## Trait.typeSurvivorship:Gradient.categoryResource ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q2.abs$beta, sei = Q2.abs$se, ci.lb = Q2.abs$ci.lb, ci.ub = Q2.abs$ci.ub,
annotate=TRUE, showweights=T, header=F,
slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
"CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
"RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))
Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts.
# Environment
# Prevalence vs Fecundity
anova(Q2.abs, X=c(-1,0,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.2751 0.1961 -1.4025 0.1608
##
## Test of Hypothesis:
## QM(df = 1) = 1.9670, p-val = 0.1608# Prevalence vs Survivorship
anova(Q2.abs, X=c(-1,0,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.0245 0.0946 -0.2592 0.7955
##
## Test of Hypothesis:
## QM(df = 1) = 0.0672, p-val = 0.7955# Intensity vs Fecundity
anova(Q2.abs, X=c(0,-1,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.4312 0.1874 -2.3008 0.0214
##
## Test of Hypothesis:
## QM(df = 1) = 5.2938, p-val = 0.0214# Intensity vs Survivorship
anova(Q2.abs, X=c(0,-1,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.1807 0.0778 -2.3232 0.0202
##
## Test of Hypothesis:
## QM(df = 1) = 5.3972, p-val = 0.0202# Pollution
# Prevalence vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,-1,0,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: 0.1805 0.2663 0.6778 0.4979
##
## Test of Hypothesis:
## QM(df = 1) = 0.4594, p-val = 0.4979# Prevalence vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,-1,0,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.1127 0.2035 -0.5540 0.5795
##
## Test of Hypothesis:
## QM(df = 1) = 0.3070, p-val = 0.5795# Intensity vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,-1,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: 0.2281 0.1902 1.1993 0.2304
##
## Test of Hypothesis:
## QM(df = 1) = 1.4383, p-val = 0.2304# Intensity vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,-1,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.0651 0.1023 -0.6363 0.5246
##
## Test of Hypothesis:
## QM(df = 1) = 0.4049, p-val = 0.5246# Resources
# Prevalence vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,-1,0,1,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.2779 0.1922 1.4456 0.1483
##
## Test of Hypothesis:
## QM(df = 1) = 2.0898, p-val = 0.1483# Prevalence vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,-1,0,0,1))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.2760 0.1994 1.3844 0.1662
##
## Test of Hypothesis:
## QM(df = 1) = 1.9165, p-val = 0.1662# Intensity vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,0,-1,1,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.0107 0.1544 0.0692 0.9448
##
## Test of Hypothesis:
## QM(df = 1) = 0.0048, p-val = 0.9448# Intensity vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,1))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.0088 0.1604 0.0548 0.9563
##
## Test of Hypothesis:
## QM(df = 1) = 0.0030, p-val = 0.9563Table3 <- data.frame(Stressor_type = rep(c("Endogenous environment", "Chemical pollution","Resource limitation"), each = 4),
Response_trait = c(rep(c("Prevalence","Intensity","Fecundity", "Survivorship"), 3)),
Overall_mean = round(Q2$beta,3),
Lower_95 = round(Q2$ci.lb,3),
Upper_95 = round(Q2$ci.ub,3),
P_value = round(Q2$pval,3),
N_effects = c(length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"]),
length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"]),
length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"]),
length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"]),
length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"]),
length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"]),
length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"]),
length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"]),
length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"]),
length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"]),
length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"]),
length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"])),
N_experiments = c(length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"]))
),
N_host_taxa = c(length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"])),
length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"]))
))
rownames(Table3) <- NULL
Table3 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | N_effects | N_experiments | N_host_taxa |
|---|---|---|---|---|---|---|---|---|
| Endogenous environment | Prevalence | 0.101 | -0.152 | 0.354 | 0.435 | 76 | 28 | 22 |
| Endogenous environment | Intensity | 0.398 | 0.214 | 0.582 | 0.000 | 152 | 55 | 41 |
| Endogenous environment | Fecundity | 0.032 | -0.517 | 0.581 | 0.910 | 9 | 4 | 4 |
| Endogenous environment | Survivorship | -0.529 | -0.705 | -0.352 | 0.000 | 162 | 66 | 47 |
| Chemical pollution | Prevalence | -0.523 | -1.017 | -0.029 | 0.038 | 19 | 9 | 7 |
| Chemical pollution | Intensity | -0.019 | -0.267 | 0.229 | 0.879 | 74 | 30 | 20 |
| Chemical pollution | Fecundity | 0.089 | -0.434 | 0.612 | 0.740 | 16 | 7 | 4 |
| Chemical pollution | Survivorship | -0.263 | -0.521 | -0.004 | 0.046 | 70 | 32 | 22 |
| Resource limitation | Prevalence | -0.150 | -0.633 | 0.333 | 0.542 | 15 | 5 | 4 |
| Resource limitation | Intensity | -0.435 | -0.792 | -0.077 | 0.017 | 35 | 13 | 10 |
| Resource limitation | Fecundity | -0.742 | -1.138 | -0.346 | 0.000 | 31 | 8 | 5 |
| Resource limitation | Survivorship | 0.005 | -0.379 | 0.388 | 0.980 | 27 | 10 | 9 |
RVE
We check if our results are qualitatively similar if we use the robust variance estimator instead of the var-covar matrix for effects estimated in the same experiment.
Q2.R <-
rma.mv(
g ~ Trait.type:Gradient.category -1,
V = dat_2$var.g,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
Q2.R_est <- clubSandwich::coef_test(Q2.R, vcov = "CR2")
forest.default(x= Q2.R_est$beta, sei = Q2.R_est$SE, pch = 19,
annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"),
slab = c("EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"),
xlab = "Standardised mean difference (Hedge's g)")
Heterogeneity
Here I take the “best” model for each question to estimate between study heterogeneity. We could potentially do the same a different model or more models, depending on our goals.
Heterogeneity is important to know if the pooled effect size can be interpreted in a meaningful way. High heterogeneity is not necessarily a bad thing, but we would at least want to know if individual study effects usually go in the direction of the pooled effect.
I2 statistic
Is one of the most commonly reported measures of heterogeneity. It captures the “percentage of variability in the effect sizes that is not caused by sampling error”. That is, heterogeneity in true effect sizes. It tells us how much true heterogeneity there is relative to sampling error.
I2 formula
#' @title Parametric simulation
#' @description Function for calculating I2 estimates using parametric simulations of model estimates taken from metafor. Note that the effectiveness of these simulations depends on the accuracy of model variance estimates.
#' @param estimate The estimate (i.e. variance) from a metafor model
#' @param sims The number of simulations
#' @param n The sample size used in estimating the variance
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @export
simMonteCarlo <- function(estimate, n, sims){
tmp <- data.frame(num = base::rep(1:sims, each = n),
y = stats::rnorm(n*sims, 0, base::sqrt(estimate)))
Var <- tmp %>% dplyr::group_by(num) %>% dplyr::summarise(Mean_var = stats::var(y))
return(as.numeric(Var$Mean_var))
}
## NOTE about PIPE: Run usethis::use_pipe() in the console. The package usethis will add what you need to import the pipe to your NAMESPACE and it will also drop warnings in checks
# Function for rounding a data frame
round_df <- function(x, digits) {
numeric_columns <- sapply(x, class) == 'numeric'
x[numeric_columns] <- round(x[numeric_columns], digits)
x
}
# Function for estimating I2
I2 <- function(model, v, ME = FALSE, sims = 1500, phylo = FALSE, obs = FALSE){
if(class(model) != "rma.mv" && class(model) != "rma"){
stop("The model object is not of class 'metafor'")
}
wi <- 1/v #weight
Vw <- sum((wi) * (length(wi) - 1)) / (((sum(wi)^2) - sum((wi)^2)))
if("rma.mv" %in% class(model) | "rma" %in% class(model)){
# Monte Carlo Simulations
# From metafor extract the important statistics
sigma2 <- matrix(model$sigma2, nrow = 1, ncol = length(model$sigma2))
colnames(sigma2) <- model$s.names
sigmaN <- model$s.nlevels
if(obs == FALSE){
stop("Please add the name of the observation-level random effect, obs. If models do not include this, re-run models including (~1|obs) in the random effect list")
}
#For each variance estimate use Monte Carlo simulation of data
Sims <- data.frame(mapply(function(x,y) simMonteCarlo(x, y, sims = sims), x = sigma2, y = sigmaN))
colnames(Sims) <- colnames(sigma2)
#Calculate total variance
VT <- rowSums(cbind(Sims, Vw))
Vt <- rowSums(Sims) # remove Vw
# For each variance component divide by the total variance. Note this needs to be fixed for phylo, but does deal with variable random effects.
I2_re <- Sims / VT
I2_total <- data.frame(Vt / VT)
tmpMatrix <- data.frame(I2_re[, -match("ID", colnames(I2_re))], total = I2_total)
names(tmpMatrix) = c(colnames(I2_re)[!colnames(I2_re) %in% 'ID'], 'total')
CI <- lapply(tmpMatrix, function(x) stats::quantile(x, c(0.025, 0.975), na.rm = TRUE))
I_CI <- as.data.frame(do.call(rbind, CI))
colnames(I_CI) <- c("2.5% CI", "97.5% CI")
I2_table <- cbind(I2_Est. = colMeans(tmpMatrix), I_CI )
class(I2_table) <- c("metaAidR", "data.frame")
return(round_df(I2_table, digits = 4))
}
}Question 1
First I calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix
I2.Q1 <- I2(Q1.R, v = dat_1$var.g, obs = "ID")*100## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'I2.Q1## I2_Est. 2.5% CI 97.5% CI
## Experiment 53.68 44.58 61.71
## total 93.86 92.71 94.86Question 2
First I calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix
I2.Q2 <- I2(Q2.R, v = dat_2$var.g, obs = "ID")*100## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'I2.Q2## I2_Est. 2.5% CI 97.5% CI
## Experiment 26.22 20.97 31.31
## total 90.27 89.31 91.13Publication bias
The most common way to assess publication bias is to test whether studies with large standard errors/low precision/small sample size tend to have larger effect sizes. Small studies have larger SEs, and because they are small, they are more likely to get published only if results are significant. In contrast, large studies, which have higher precision and have costed more money, are more likely to get published either way. For this reason small studies can cause publication bias.
Funnel plots
A very common way to assess publication bias is therefore to check whether studies with low precision generally have larger effects. The typical way to visualize this is through a funnel plot with all effects. Now, this is not entirely correct for us, as the effects are not independent. Nakagawa et al 2021 (Methods Ecol Evol) suggest using conditional residuals for funnel plots. These residual plots take some of the non-independence into account, but are not perfect, as they assume: 1)sampling SE(sei) does not covary with moderators in meta-regression, 2) sampling SE is the same as the SE of the residuals and 3) sampling variances are not correlated. This last one is most certainly incorrect.
Estimation of conditional residuals is not available to rma.mv objects (which can include a sampling variance-covariance matrix). The available alternative is to ignore this covariation in sampling variance for the assessment of publication bias. This is a caveat that should be acknowledgend in the paper, yet I think it is a useful visualization nonetheless. To do these funnel plots I re-fit the best model as rma.uni for each data set. While not ideal, this the best current alternative suggested by Nakagawa et al 2021 (Methods Ecol Evol)
Note I checked some previous studies that use multilevel regression models as we do here, to see what they do about publication bias. Sauer et al 2020 Ecology (101) acknowledged this problem and decided not to even show funnel plots: “Because of the complex non independence among effect sizes within a study (e.g., some studies had multiple effect sizes), we did not use funnel plots or rank correlation tests to assess publication bias (Lau et al. 2006, Civitello et al.2015).” Others (Yates et al. 2019 Environmental DNA (1); Salerno et al 2021 Global Change Biology (27); Brlik et al. 2020 JAE) acknowledge this problem, but still show funnel plots warning caution. They also conduct a version of Eggerts test somewhat in the same spirit as I do below.
Data set 1: fitness effects of environmental stress
# re fit as a mixed effects model with rma.uni
# any fitting method other than "FE" works for a random/mixed effects model
Q1.uni <-
rma.uni(
g ~ Trait.type:Gradient.category + Experiment + ID ,
vi = dat_1$var.g,
data = dat_1,
method = "REML"
)## Warning: Redundant predictors dropped from the model.# get conditional residuals (residuals for each effect)
Q1res <- rstandard.rma.uni(Q1.uni, type = "conditional")
# create colour palette as gradient of effect size g
nlvl <- length(dat_1$g)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))(nlvl)
pal <-pal[order(match(pal,dat_1$g))]
# plot funnel with conditional residuals
metafor::funnel(Q1res$resid, sqrt(dat_1$var.g), level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), col=pal, back = "white", refline=0, legend=TRUE, xlab = "Conditional residual", ylab = "Standard error")
We can see that as SE increases (e.g ~ 1) all studies report large negative fitness effects of environmental stress (large negative residuals). In other words, it seems that very small studies tend to find that environmental stress is more negative to fitness!
Data set 2: Fitness and infectivity effects of environmental stress
# re fit as a mixed effects model with rma.uni
# any fitting method other than "FE" works for a random/mixed effects model
Q2.uni <-
rma.uni(
g ~ Trait.type:Gradient.category + Experiment + ID,
vi = dat_2$var.g,
data = dat_2,
method = "REML"
)## Warning: Redundant predictors dropped from the model.# get conditional residuals (residuals for each effect)
Q2res <- rstandard.rma.uni(Q2.uni, type = "conditional")
# create colour palette as gradient of effect size g
nlvl <- length(dat_2$g)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))(nlvl)
pal <-pal[order(match(pal,dat_2$g))]
# plot funnel with all conditional residuals
metafor::funnel(Q2res$resid, sqrt(dat_2$var.g), level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), col=pal, back = "white", refline=0, legend=TRUE, xlab = "Conditional residual", ylab = "Standard error")
Because this data set includes both the infectivity and fitness effects it’s not surprising that small studies are distributed on both sides of the funnel. This does not mean absence of publication bias as negative effects for fitness and positive effects for (-) infectivity both would suggest that small studies tend to inflate the negative consequences of environmental stress. For the publication, I suggest we use -g for infectivity effects (currently g = infection intensity or prevalence) and/or we use funnel plots for only dataset 1.
Two step approach to test for publication bias
Following Nakagawa et al. (2021) we conduct a two-step approach to test for publication bias, that can be thought of as a modified Egger’s test for multilevel meta analysis. First we include the effect sizes’ standard error as the only moderator and including the same random effect structure as in our meta-analysis. If the slope of this moderator is significant, the estimated intercept provides a downwardly biased estimate of the overall meta-analytic effect (Doucouliagos 2012, 2014). In this case, we run a second meta regression with only the effect sizes’ variances as moderators. The intercept in this case may still be biased, but less so.
Data set 1: fitness effects of environmental stress
# get standard errors of effect sizes
dat_1$sei <- sqrt(dat_1$var.g)
# meta regression with SE
d1s1 <- rma.mv(
g ~ 1 + sei,
V = varcovmat_1_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
d1s1##
## Multivariate Meta-Analysis Model (k = 384; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4132 0.6428 78 no Experiment
## sigma^2.2 0.2912 0.5396 384 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 382) = 247179.2283, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 56.3837, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3606 0.1412 2.5528 0.0107 0.0837 0.6374 *
## sei -1.5464 0.2059 -7.5089 <.0001 -1.9501 -1.1428 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The slope is negative and significant. Small studies report large and negative effect sizes!! Because the slope of sei is significant, variance is a better moderator according to (Doucouliagos 2012, 2014).
d1s2 <- rma.mv(
g ~ 1 + dat_1$var.g,
V = varcovmat_1_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
d1s2 ##
## Multivariate Meta-Analysis Model (k = 384; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3473 0.5893 78 no Experiment
## sigma^2.2 0.2943 0.5425 384 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 382) = 422619.1025, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 77.9257, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.1613 0.0911 -1.7708 0.0766 -0.3398 0.0172 .
## dat_1$var.g -0.8711 0.0987 -8.8276 <.0001 -1.0645 -0.6777 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Again, precision effects are significant, suggesting publication bias. Thus, to “account” for potential effects of publication bias, we need to include variance as a moderator.
Including the variance-covariance error matrix within studies
Q1.B <-
rma.mv(
g ~ Trait.type:Gradient.category + var.g -1,
V = varcovmat_1_PD$mat,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)
summary(Q1.B)##
## Multivariate Meta-Analysis Model (k = 384; method: REML)
##
## logLik Deviance AIC BIC AICc
## -507.8909 1015.7818 1033.7818 1069.1720 1034.2723
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3149 0.5612 78 no Experiment
## sigma^2.2 0.2826 0.5316 384 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 377) = 133298.1893, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 125.9072, p-val < .0001
##
## Model Results:
##
## estimate se zval
## var.g -0.8846 0.0980 -9.0217
## Trait.typeFecundity:Gradient.categoryEnvironment -0.0978 0.2332 -0.4196
## Trait.typeSurvivorship:Gradient.categoryEnvironment -0.1762 0.1109 -1.5890
## Trait.typeFecundity:Gradient.categoryPollution 0.0589 0.3121 0.1886
## Trait.typeSurvivorship:Gradient.categoryPollution -0.1061 0.1548 -0.6853
## Trait.typeFecundity:Gradient.categoryResource -0.6939 0.2177 -3.1877
## Trait.typeSurvivorship:Gradient.categoryResource 0.3048 0.2032 1.5002
## pval ci.lb ci.ub
## var.g <.0001 -1.0767 -0.6924
## Trait.typeFecundity:Gradient.categoryEnvironment 0.6747 -0.5548 0.3591
## Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1121 -0.3934 0.0411
## Trait.typeFecundity:Gradient.categoryPollution 0.8504 -0.5529 0.6706
## Trait.typeSurvivorship:Gradient.categoryPollution 0.4932 -0.4094 0.1973
## Trait.typeFecundity:Gradient.categoryResource 0.0014 -1.1206 -0.2673
## Trait.typeSurvivorship:Gradient.categoryResource 0.1336 -0.0934 0.7030
##
## var.g ***
## Trait.typeFecundity:Gradient.categoryEnvironment
## Trait.typeSurvivorship:Gradient.categoryEnvironment
## Trait.typeFecundity:Gradient.categoryPollution
## Trait.typeSurvivorship:Gradient.categoryPollution
## Trait.typeFecundity:Gradient.categoryResource **
## Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q1.B$beta, sei = Q1.B$se, pch = 19,
annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"),
slab = c("variance", "EE:fecundity", "EE:survivorship", "CP:fecundity", "CP:survivorship", "RL:fecundity", "RL:survivorship"),
xlab = "Standardised mean difference (Hedge's g)")
Data set 2: fitness and infectivity effects of environmental stress
# get standard errors of effect sizes
dat_2$sei <- sqrt(dat_2$var.g)
# meta regression with SE
d2s1 <- rma.mv(
g ~ 1 + sei,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML"
)As in data set 1, effects of small studies are dominated by negative effect sizes. Because the slope of sei is significant the variance is a better moderator for the overall effect.
# meta regression with SE
d2s2 <- rma.mv(
g ~ 1 + dat_2$var.g,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML"
)
d2s2##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1407 0.3751 113 no Experiment
## sigma^2.2 0.5141 0.7170 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 684) = 3689465.0418, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 16.9025, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.0081 0.0597 -0.1363 0.8916 -0.1251 0.1088
## dat_2$var.g -0.3665 0.0891 -4.1113 <.0001 -0.5411 -0.1918 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Q2.B <-
rma.mv(
g ~ Trait.type:Gradient.category + var.g -1,
V = varcovmat_2_PD$mat,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML"
)
summary(Q2.B)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -1020.1000 2040.2001 2070.2001 2137.8762 2070.9307
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1691 0.4112 113 no Experiment
## sigma^2.2 0.4257 0.6525 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 673) = 827885.4316, p-val < .0001
##
## Test of Moderators (coefficients 1:13):
## QM(df = 13) = 124.3757, p-val < .0001
##
## Model Results:
##
## estimate se zval
## var.g -0.3826 0.0884 -4.3297
## Trait.typePrevalence:Gradient.categoryEnvironment 0.2548 0.1326 1.9217
## Trait.typeIntensity:Gradient.categoryEnvironment 0.5017 0.0958 5.2366
## Trait.typeFecundity:Gradient.categoryEnvironment 0.1025 0.2754 0.3723
## Trait.typeSurvivorship:Gradient.categoryEnvironment -0.4050 0.0932 -4.3449
## Trait.typePrevalence:Gradient.categoryPollution -0.3950 0.2506 -1.5764
## Trait.typeIntensity:Gradient.categoryPollution 0.0501 0.1261 0.3970
## Trait.typeFecundity:Gradient.categoryPollution 0.2032 0.2645 0.7685
## Trait.typeSurvivorship:Gradient.categoryPollution -0.1509 0.1326 -1.1378
## Trait.typePrevalence:Gradient.categoryResource -0.0399 0.2440 -0.1635
## Trait.typeIntensity:Gradient.categoryResource -0.3309 0.1816 -1.8224
## Trait.typeFecundity:Gradient.categoryResource -0.6764 0.1999 -3.3847
## Trait.typeSurvivorship:Gradient.categoryResource 0.1080 0.1942 0.5562
## pval ci.lb ci.ub
## var.g <.0001 -0.5557 -0.2094
## Trait.typePrevalence:Gradient.categoryEnvironment 0.0546 -0.0051 0.5146
## Trait.typeIntensity:Gradient.categoryEnvironment <.0001 0.3139 0.6895
## Trait.typeFecundity:Gradient.categoryEnvironment 0.7097 -0.4373 0.6424
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 -0.5877 -0.2223
## Trait.typePrevalence:Gradient.categoryPollution 0.1149 -0.8862 0.0961
## Trait.typeIntensity:Gradient.categoryPollution 0.6913 -0.1971 0.2973
## Trait.typeFecundity:Gradient.categoryPollution 0.4422 -0.3151 0.7216
## Trait.typeSurvivorship:Gradient.categoryPollution 0.2552 -0.4108 0.1090
## Trait.typePrevalence:Gradient.categoryResource 0.8701 -0.5182 0.4384
## Trait.typeIntensity:Gradient.categoryResource 0.0684 -0.6868 0.0250
## Trait.typeFecundity:Gradient.categoryResource 0.0007 -1.0681 -0.2847
## Trait.typeSurvivorship:Gradient.categoryResource 0.5781 -0.2726 0.4886
##
## var.g ***
## Trait.typePrevalence:Gradient.categoryEnvironment .
## Trait.typeIntensity:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryEnvironment
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typePrevalence:Gradient.categoryPollution
## Trait.typeIntensity:Gradient.categoryPollution
## Trait.typeFecundity:Gradient.categoryPollution
## Trait.typeSurvivorship:Gradient.categoryPollution
## Trait.typePrevalence:Gradient.categoryResource
## Trait.typeIntensity:Gradient.categoryResource .
## Trait.typeFecundity:Gradient.categoryResource ***
## Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q2.B$beta, sei = Q2.B$se, pch = 19,
annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"),
slab = c("variance","EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"),
xlab = "Standardised mean difference (Hedge's g)")
Contrasts
Re-run model with absolute values
Q2.B.abs <-
rma.mv(
absolute ~ Trait.type:Gradient.category + var.g -1,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
summary(Q2.B.abs)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -671.1227 1342.2454 1372.2454 1439.9216 1372.9760
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0890 0.2983 113 no Experiment
## sigma^2.2 0.0712 0.2669 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 673) = 822493.8021, p-val < .0001
##
## Test of Moderators (coefficients 1:13):
## QM(df = 13) = 633.2958, p-val < .0001
##
## Model Results:
##
## estimate se zval
## var.g 1.0912 0.0709 15.4002
## Trait.typePrevalence:Gradient.categoryEnvironment 0.4348 0.0888 4.8977
## Trait.typeIntensity:Gradient.categoryEnvironment 0.6922 0.0661 10.4710
## Trait.typeFecundity:Gradient.categoryEnvironment 0.3123 0.1641 1.9036
## Trait.typeSurvivorship:Gradient.categoryEnvironment 0.5140 0.0625 8.2278
## Trait.typePrevalence:Gradient.categoryPollution 0.3522 0.1796 1.9608
## Trait.typeIntensity:Gradient.categoryPollution 0.4506 0.0850 5.2989
## Trait.typeFecundity:Gradient.categoryPollution 0.5483 0.1760 3.1149
## Trait.typeSurvivorship:Gradient.categoryPollution 0.3254 0.0879 3.7015
## Trait.typePrevalence:Gradient.categoryResource 0.2295 0.1593 1.4409
## Trait.typeIntensity:Gradient.categoryResource 0.4900 0.1185 4.1343
## Trait.typeFecundity:Gradient.categoryResource 0.5740 0.1297 4.4262
## Trait.typeSurvivorship:Gradient.categoryResource 0.5195 0.1231 4.2215
## pval ci.lb ci.ub
## var.g <.0001 0.9524 1.2301
## Trait.typePrevalence:Gradient.categoryEnvironment <.0001 0.2608 0.6088
## Trait.typeIntensity:Gradient.categoryEnvironment <.0001 0.5626 0.8218
## Trait.typeFecundity:Gradient.categoryEnvironment 0.0570 -0.0093 0.6338
## Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001 0.3916 0.6365
## Trait.typePrevalence:Gradient.categoryPollution 0.0499 0.0002 0.7041
## Trait.typeIntensity:Gradient.categoryPollution <.0001 0.2839 0.6172
## Trait.typeFecundity:Gradient.categoryPollution 0.0018 0.2033 0.8933
## Trait.typeSurvivorship:Gradient.categoryPollution 0.0002 0.1531 0.4977
## Trait.typePrevalence:Gradient.categoryResource 0.1496 -0.0827 0.5417
## Trait.typeIntensity:Gradient.categoryResource <.0001 0.2577 0.7223
## Trait.typeFecundity:Gradient.categoryResource <.0001 0.3198 0.8281
## Trait.typeSurvivorship:Gradient.categoryResource <.0001 0.2783 0.7607
##
## var.g ***
## Trait.typePrevalence:Gradient.categoryEnvironment ***
## Trait.typeIntensity:Gradient.categoryEnvironment ***
## Trait.typeFecundity:Gradient.categoryEnvironment .
## Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## Trait.typePrevalence:Gradient.categoryPollution *
## Trait.typeIntensity:Gradient.categoryPollution ***
## Trait.typeFecundity:Gradient.categoryPollution **
## Trait.typeSurvivorship:Gradient.categoryPollution ***
## Trait.typePrevalence:Gradient.categoryResource
## Trait.typeIntensity:Gradient.categoryResource ***
## Trait.typeFecundity:Gradient.categoryResource ***
## Trait.typeSurvivorship:Gradient.categoryResource ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q2.B.abs$beta, sei = Q2.B.abs$se, pch = 19,
annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"),
slab = c("variance","EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"),
xlab = "Standardised mean difference (Hedge's g)")
Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts.
# Environment
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,-1,0,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.1225 0.1768 -0.6930 0.4883
##
## Test of Hypothesis:
## QM(df = 1) = 0.4802, p-val = 0.4883# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,-1,0,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: 0.0792 0.0865 0.9156 0.3599
##
## Test of Hypothesis:
## QM(df = 1) = 0.8383, p-val = 0.3599# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,-1,1,0,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.3799 0.1684 -2.2559 0.0241
##
## Test of Hypothesis:
## QM(df = 1) = 5.0889, p-val = 0.0241# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,-1,0,1,0,0,0,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0
##
## Results:
## estimate se zval pval
## 1: -0.1782 0.0716 -2.4885 0.0128
##
## Test of Hypothesis:
## QM(df = 1) = 6.1928, p-val = 0.0128# Pollution
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,-1,0,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: 0.1962 0.2431 0.8070 0.4197
##
## Test of Hypothesis:
## QM(df = 1) = 0.6512, p-val = 0.4197# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,-1,0,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.0267 0.1856 -0.1441 0.8854
##
## Test of Hypothesis:
## QM(df = 1) = 0.0208, p-val = 0.8854# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,-1,1,0,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: 0.0977 0.1747 0.5593 0.5759
##
## Test of Hypothesis:
## QM(df = 1) = 0.3129, p-val = 0.5759# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,-1,0,1,0,0,0,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0
##
## Results:
## estimate se zval pval
## 1: -0.1252 0.0937 -1.3358 0.1816
##
## Test of Hypothesis:
## QM(df = 1) = 1.7843, p-val = 0.1816# Resources
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,1,0))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.3445 0.1754 1.9640 0.0495
##
## Test of Hypothesis:
## QM(df = 1) = 3.8574, p-val = 0.0495# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,0,1))##
## Hypothesis:
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.2900 0.1826 1.5885 0.1122
##
## Test of Hypothesis:
## QM(df = 1) = 2.5233, p-val = 0.1122# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,0,-1,1,0))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.0840 0.1405 0.5975 0.5502
##
## Test of Hypothesis:
## QM(df = 1) = 0.3570, p-val = 0.5502# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,0,-1,0,1))##
## Hypothesis:
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0
##
## Results:
## estimate se zval pval
## 1: 0.0295 0.1452 0.2031 0.8391
##
## Test of Hypothesis:
## QM(df = 1) = 0.0412, p-val = 0.8391Supporting Tables
Table S3
TableS3 <- data.frame(question = rep(c("Q1", "Q2"), each = 2),
step = rep(c("Step 1", "Step 2"),2),
slope = c(round(d1s1$b[2],3), round(d1s2$b[2],3), round(d2s1$b[2],3), round(d2s2$b[2],3)),
se = c(round(d1s1$se[2],3), round(d1s2$se[2],3), round(d2s1$se[2],3), round(d2s2$se[2],3)),
pval =c(round(d1s1$pval[2],3), round(d1s2$pval[2],3), round(d2s1$pval[2],3), round(d2s2$pval[2],3)))
TableS3 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| question | step | slope | se | pval |
|---|---|---|---|---|
| Q1 | Step 1 | -1.546 | 0.206 | 0.000 |
| Q1 | Step 2 | -0.871 | 0.099 | 0.000 |
| Q2 | Step 1 | -0.544 | 0.165 | 0.001 |
| Q2 | Step 2 | -0.366 | 0.089 | 0.000 |
Table S4
TableS4 <- data.frame(Stressor_type = Table2$Stressor_type,
Response_trait = Table2$Response_trait,
Overall_mean = Table2$Overall_mean,
Lower_95 = Table2$Lower_95,
Upper_95 = Table2$Upper_95,
P_value = Table2$P_value,
Overall_mean_adjusted = round(Q1.B$beta[-1],3),
Lower_95_adjusted = round(Q1.B$ci.lb[-1],3),
Upper_95_adjusted = round(Q1.B$ci.ub[-1],3),
P_value_adjusted = round(Q1.B$pval[-1],3),
N_effects = Table2$N_effects,
N_experiments = Table2$N_experiments,
N_host_taxa = Table2$N_host_taxa)
TableS4 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = T)| Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | Overall_mean_adjusted | Lower_95_adjusted | Upper_95_adjusted | P_value_adjusted | N_effects | N_experiments | N_host_taxa |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Endogenous environment | Fecundity | -0.326 | -0.802 | 0.150 | 0.180 | -0.098 | -0.555 | 0.359 | 0.675 | 14 | 3 | 3 |
| Endogenous environment | Survivorship | -0.516 | -0.715 | -0.317 | 0.000 | -0.176 | -0.393 | 0.041 | 0.112 | 189 | 45 | 34 |
| Chemical pollution | Fecundity | -0.141 | -0.759 | 0.477 | 0.656 | 0.059 | -0.553 | 0.671 | 0.850 | 14 | 4 | 2 |
| Chemical pollution | Survivorship | -0.351 | -0.640 | -0.062 | 0.017 | -0.106 | -0.409 | 0.197 | 0.493 | 76 | 23 | 15 |
| Resource limitation | Fecundity | -0.864 | -1.287 | -0.440 | 0.000 | -0.694 | -1.121 | -0.267 | 0.001 | 59 | 6 | 4 |
| Resource limitation | Survivorship | 0.002 | -0.399 | 0.404 | 0.990 | 0.305 | -0.093 | 0.703 | 0.134 | 32 | 6 | 5 |
Table S5
TableS5 <- data.frame(Stressor_type = Table3$Stressor_type,
Response_trait = Table3$Response_trait,
Overall_mean = Table3$Overall_mean,
Lower_95 = Table3$Lower_95,
Upper_95 = Table3$Upper_95,
P_value = Table3$P_value,
Overall_mean_adjusted = round(Q2.B$beta[-1],3),
Lower_95_adjusted = round(Q2.B$ci.lb[-1],3),
Upper_95_adjusted = round(Q2.B$ci.ub[-1],3),
P_value_adjusted = round(Q2.B$pval[-1],3),
N_effects = Table3$N_effects,
N_experiments = Table3$N_experiments,
N_host_taxa = Table3$N_host_taxa)
TableS5 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = T)| Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | Overall_mean_adjusted | Lower_95_adjusted | Upper_95_adjusted | P_value_adjusted | N_effects | N_experiments | N_host_taxa |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Endogenous environment | Prevalence | 0.101 | -0.152 | 0.354 | 0.435 | 0.255 | -0.005 | 0.515 | 0.055 | 76 | 28 | 22 |
| Endogenous environment | Intensity | 0.398 | 0.214 | 0.582 | 0.000 | 0.502 | 0.314 | 0.689 | 0.000 | 152 | 55 | 41 |
| Endogenous environment | Fecundity | 0.032 | -0.517 | 0.581 | 0.910 | 0.103 | -0.437 | 0.642 | 0.710 | 9 | 4 | 4 |
| Endogenous environment | Survivorship | -0.529 | -0.705 | -0.352 | 0.000 | -0.405 | -0.588 | -0.222 | 0.000 | 162 | 66 | 47 |
| Chemical pollution | Prevalence | -0.523 | -1.017 | -0.029 | 0.038 | -0.395 | -0.886 | 0.096 | 0.115 | 19 | 9 | 7 |
| Chemical pollution | Intensity | -0.019 | -0.267 | 0.229 | 0.879 | 0.050 | -0.197 | 0.297 | 0.691 | 74 | 30 | 20 |
| Chemical pollution | Fecundity | 0.089 | -0.434 | 0.612 | 0.740 | 0.203 | -0.315 | 0.722 | 0.442 | 16 | 7 | 4 |
| Chemical pollution | Survivorship | -0.263 | -0.521 | -0.004 | 0.046 | -0.151 | -0.411 | 0.109 | 0.255 | 70 | 32 | 22 |
| Resource limitation | Prevalence | -0.150 | -0.633 | 0.333 | 0.542 | -0.040 | -0.518 | 0.438 | 0.870 | 15 | 5 | 4 |
| Resource limitation | Intensity | -0.435 | -0.792 | -0.077 | 0.017 | -0.331 | -0.687 | 0.025 | 0.068 | 35 | 13 | 10 |
| Resource limitation | Fecundity | -0.742 | -1.138 | -0.346 | 0.000 | -0.676 | -1.068 | -0.285 | 0.001 | 31 | 8 | 5 |
| Resource limitation | Survivorship | 0.005 | -0.379 | 0.388 | 0.980 | 0.108 | -0.273 | 0.489 | 0.578 | 27 | 10 | 9 |
R1
Here we conduct analyses suggested by reviewers.
Correct for sample size and get g
# sample size correction factor
for (i in 1:N) {
# if odds ratio or comparison between means
if (is.na(dat_UI$N_U[i]) == F) {
dat_UI$J[i] <- 1 - 3 / (4 * (dat_UI$N_U[i] + dat_UI$N_I[i] - 2) - 1)
} else {
# if correlation
if (is.na(dat_UI$R_N[i]) == F) {
dat_UI$J[i] <- 1 - 3 / (4 * (dat_UI$R_N[i] - 2) - 1)
}else {
dat_UI$J[i] <- NA
}
}
}
# corrected effect size
dat_UI$g <- dat_UI$J * dat_UI$d
# and variance
dat_UI$var.g <- (dat_UI$J ^ 2) * dat_UI$var.d
# remove studies without variance for which g had been set to NA
dat_UI <- dat_UI[!is.na(dat_UI$g),]
# remove duplicated rows which were created whenever controls were used more than once in main analysis.
dat_UI <- dat_UI[!duplicated(dat_UI[,-c(1,21)]), ]Flip effect sizes for remaining mortality effects
mort <- grep(pattern = "[M,m]ort", x = dat_UI$Trait.x)
for(i in mort){
dat_UI$g[i] <- -dat_UI$g[i]
}Identify negative eigenvalues
As in our main analysis, we find out which studies have negative eigenvalues in their variance-covar matrix. This means covariance between treatment effects with the same control is larger than variance and it’s an indicator of suspiciously low variance and potentially an incorrectly reported N or an error while digitizing data.
studies_1 <- unique(dat_UI$Study)
studies_to_Uheck <- c()
negative_eigen <- c()
for (k in studies_1) {
dat_study <- dat_UI[which(dat_UI$Study == k),]
varcovmat = matrix(0,
nrow = dim(dat_study)[1],
ncol = dim(dat_study)[1])
for (i in 1:dim(dat_study)[1]) {
for (j in 1:dim(dat_study)[1]) {
if (i == j) {
varcovmat[i, j] = dat_study$var.g[i]
} else{
if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
dat_study[i, "Level_U"] == dat_study[j, "Level_U"]) {
varcovmat[i, j] = 1 / dat_study[i, "N_U"] + dat_study$g[i] * dat_study$g[j] /
(dat_study[i, "N_U"] + dat_study[i, "N_I"] + dat_study[j, "N_I"])
}
}
}
}
val <- eigen(varcovmat)
for (m in 1:length(val$values)) {
if (val$values[m] < 0) {
studies_to_Uheck <- append(studies_to_Uheck, k)
negative_eigen <- append(negative_eigen, val$values[m])
}
}
}
Neg_eigen_UI <- data.frame(study = studies_to_Uheck, eigen = negative_eigen)
Neg_eigen_UI %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| study | eigen |
|---|---|
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.0181032 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.1759683 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.3098886 |
| Ashraf et al. 2017 Environ Sci Pollut Res | -0.8423039 |
| Manzi et al. 2019. Freshwater Biology | -0.0103166 |
| Manzi et al. 2019. Freshwater Biology | -0.0327272 |
# exclude suspicious studies with negative eigenvalues and which already raised flags in the main analysis
dat_UI <- dat_UI[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_UI$Study),]
dat_UI <- dat_UI[-grep("Civitello et al. 2020. Proc B", dat_UI$Study),]
# I'm also excluding Hall et al. as they also had too negative eigen values
dat_UI <- dat_UI[-grep("Hall et al. 2013. Oecologia", dat_UI$Study),]
dat_UI <- dat_UI[-grep("Manzi et al. 2019. Freshwater Biology", dat_UI$Study),]
varcovmat_UI = matrix(0, nrow = dim(dat_UI)[1], ncol = dim(dat_UI)[1])
for (i in 1:dim(dat_UI)[1]) {
for (j in 1:dim(dat_UI)[1]) {
if (i == j) {varcovmat_UI[i,j] = dat_UI$var.g[i]}else{
if (dat_UI[i, "Experiment"] == dat_UI[j, "Experiment"] & dat_UI[i, "Level_U"] == dat_UI[j, "Level_U"]) {
varcovmat_UI[i,j] = 1/dat_UI[i,"N_U"] + dat_UI$g[i]*dat_UI$g[j]/(dat_UI[i,"N_U"] + dat_UI[i,"N_I"] + dat_UI[j,"N_I"])
}
}
}
}
# correct negative eigens in the few studies with large covar to var ratios
varcovmat_UI_PD <- nearPD(varcovmat_UI, keepDiag = TRUE)Model without moderators
We fit the model to test whether infection reduces fitness (fecundity or survivorship)
Q0.R1.f <-
rma.mv(
g ~ Trait.type -1,
V = varcovmat_UI_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID.x),
data = dat_UI,
method = "REML"
)
pal <- RColorBrewer::brewer.pal(6,"Spectral")
Q0.R1.f_OP <-
orchard_plot(
Q0.R1.f,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
scale_fill_manual(values = rep(pal, each = 2)) + scale_colour_manual(values = rep(pal, each = 2)) +
theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_y_discrete(
labels = c(
"Fecundity",
"Survivorship"
)
)
Q0.R1.f_OP
Forest plot
summary(Q0.R1.f)##
## Multivariate Meta-Analysis Model (k = 308; method: REML)
##
## logLik Deviance AIC BIC AICc
## -1161.8670 2323.7340 2331.7340 2346.6284 2331.8669
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.9933 0.9966 72 no Experiment
## sigma^2.2 0.0259 0.1608 118 no ID.x
##
## Test for Residual Heterogeneity:
## QE(df = 306) = 3629.6611, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 56.5456, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## Trait.typeFecundity -2.0056 0.2984 -6.7215 <.0001 -2.5905 -1.4208
## Trait.typeSurvivorship -0.6775 0.1428 -4.7457 <.0001 -0.9573 -0.3977
##
## Trait.typeFecundity ***
## Trait.typeSurvivorship ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(x= Q0.R1.f$beta, sei = Q0.R1.f$se, ci.lb = Q0.R1.f$ci.lb, ci.ub = Q0.R1.f$ci.ub,
annotate=TRUE, showweights=T, header=F,
slab = c("Fecundity", "Survivorship"))
Comment 2.R1 - Do aquatic and terrestrial systems differ in their responses to stressors?
To address this question by Reviewer 2, we’ll add an ‘Environment’ moderator to the final models in Q1 and Q2, presented in the main text. The results of these models are presented in the Supporting Material.
Q1: Best Model + Environment
Effects in survival and fecundity under different types of stress and in aquatic vs terrestrial habitats, combining infected and uninfected.
Q1.R1 <-
rma.mv(
g ~ Environment:Trait.type:Gradient.category -1 ,
V = varcovmat_1_PD$mat,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)## Warning: Redundant predictors dropped from the model.summary(Q1.R1)##
## Multivariate Meta-Analysis Model (k = 384; method: REML)
##
## logLik Deviance AIC BIC AICc
## -540.8105 1081.6210 1103.6210 1146.8172 1104.3483
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2203 0.4693 78 no Experiment
## sigma^2.2 0.3630 0.6025 384 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 375) = 214780.6070, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 57.7444, p-val < .0001
##
## Model Results:
##
## estimate
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment -0.3890
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.6141
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1304
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.0724
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -0.1986
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -0.4953
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -0.9815
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -0.4179
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource -0.0528
## se
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.2419
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1048
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.2793
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.3140
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.1982
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 0.2038
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource 0.2331
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 0.4802
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.2039
## zval
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment -1.6084
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -5.8582
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.4668
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.2306
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -1.0019
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -2.4305
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -4.2102
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -0.8702
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource -0.2588
## pval
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.1078
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.6407
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.8176
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.3164
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 0.0151
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource <.0001
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 0.3842
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.7958
## ci.lb
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment -0.8630
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.8196
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.4170
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.6879
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -0.5871
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -0.8947
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -1.4384
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -1.3591
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource -0.4523
## ci.ub
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.0850
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.4086
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.6777
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.5430
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.1899
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -0.0959
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -0.5246
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 0.5233
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.3468
##
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution *
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource ***
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(
x = Q1.R1$beta,
sei = Q1.R1$se,
ci.lb = Q1.R1$ci.lb,
ci.ub = Q1.R1$ci.ub,
annotate = TRUE,
showweights = T,
header = F,
slab = c(
"Aquatic environment fecundity",
"Aquatic environment survivorship",
"Terrestrial environment survivorship",
"Aquatic pollution fecundity",
"Aquatic pollution survivorship",
"Terrestrial pollution survivorship",
"Aquatic resource fecundity",
"Terrestrial resource fecundity",
"Aquatic resource survivorship"
)
)
pal <- RColorBrewer::brewer.pal(11,"Spectral")
Q1.R1_OP <-
orchard_plot(
Q1.R1,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
scale_fill_manual(values = pal[c(7,8,8,10,11,11,5,5,4)]) +
scale_colour_manual(values = pal[c(7,8,8,10,11,11,5,5,4)]) +
theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_y_discrete(
labels = c("AQ:EE:fecundity",
"AQ:EE:survivorship",
"TE:EE:survivorship",
"AQ:CP:fecundity",
"AQ:CP:survivorship",
"TE:CP:survivorship",
"AQ:RL:fecundity",
"TE:RL:fecundity",
"AQ:RL:survivorship"))## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.Q1.R1_OP
Make Table with results from Q1.R1
TableQ1.R1 <-
data.frame(
Environment = c(
"Aquatic",
"Aquatic",
"Terrestrial",
"Aquatic",
"Aquatic",
"Terrestrial",
"Aquatic",
"Terrestrial",
"Aquatic"
),
Gradient.category = rep(
c(
"Environment",
"Pollution",
"Resource"
),
each = 3
),
Trait.type = c(
"Fecundity",
"Survivorship",
"Survivorship",
"Fecundity",
"Survivorship",
"Survivorship",
"Fecundity",
"Fecundity",
"Survivorship"
),
Overall_mean = round(Q1.R1$beta, 3),
Lower_95 = round(Q1.R1$ci.lb, 3),
Upper_95 = round(Q1.R1$ci.ub, 3),
P_value = round(Q1.R1$pval, 3))
N_effects <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_effects = length(ID))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.N_experiments <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_experiments = length(unique(Experiment)))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.N_hosts <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_hosts = length(unique(Host)))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.TableQ1.R1 <-
left_join(TableQ1.R1, N_effects) %>% left_join(N_experiments) %>% left_join(N_experiments) %>%
dplyr::rename(Stressor_type = Gradient.category,
Response_trait = Trait.type) %>%
mutate(Stressor_type = plyr::revalue(
Stressor_type,
c(
"Environment" = "Endogenous Environment",
"Pollution" = "Chemical Pollution",
"Resource" = "Resource Limitation"
)
))## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type,
## N_experiments)`rownames(TableQ1.R1) <- NULL
TableQ1.R1 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Environment | Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | N_effects | N_experiments |
|---|---|---|---|---|---|---|---|---|
| Aquatic | Endogenous Environment | Fecundity | -0.389 | -0.863 | 0.085 | 0.108 | 14 | 3 |
| Aquatic | Endogenous Environment | Survivorship | -0.614 | -0.820 | -0.409 | 0.000 | 170 | 40 |
| Terrestrial | Endogenous Environment | Survivorship | 0.130 | -0.417 | 0.678 | 0.641 | 19 | 5 |
| Aquatic | Chemical Pollution | Fecundity | -0.072 | -0.688 | 0.543 | 0.818 | 14 | 4 |
| Aquatic | Chemical Pollution | Survivorship | -0.199 | -0.587 | 0.190 | 0.316 | 44 | 12 |
| Terrestrial | Chemical Pollution | Survivorship | -0.495 | -0.895 | -0.096 | 0.015 | 32 | 11 |
| Aquatic | Resource Limitation | Fecundity | -0.981 | -1.438 | -0.525 | 0.000 | 23 | 5 |
| Terrestrial | Resource Limitation | Fecundity | -0.418 | -1.359 | 0.523 | 0.384 | 36 | 1 |
| Aquatic | Resource Limitation | Survivorship | -0.053 | -0.452 | 0.347 | 0.796 | 32 | 6 |
Does accounting for the transmission environment reduce heterogeneity? We calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix.
Q1.R1.R <-
rma.mv(
g ~ Environment:Trait.type:Gradient.category -1 ,
V = dat_1$var.g,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_1,
method = "REML"
)## Warning: Redundant predictors dropped from the model.I2.Q1.R1 <- I2(Q1.R1.R, v = dat_1$var.g, obs = "ID")*100## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'I2.Q1.R1## I2_Est. 2.5% CI 97.5% CI
## Experiment 51.52 42.40 59.52
## total 93.63 92.42 94.63Q2: Best Model + Environment
Trait type x Gradient category x Transmission environment interaction
Q2.R1 <-
rma.mv(
g ~ Environment:Trait.type:Gradient.category -1,
V = varcovmat_2_PD$mat,
random = list(~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML",
)
summary(Q2.R1)##
## Multivariate Meta-Analysis Model (k = 686; method: REML)
##
## logLik Deviance AIC BIC AICc
## -999.5082 1999.0163 2051.0163 2167.8932 2053.2273
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1790 0.4231 113 no Experiment
## sigma^2.2 0.4244 0.6515 686 no ID
##
## Test for Residual Heterogeneity:
## QE(df = 662) = 1419271.1848, p-val < .0001
##
## Test of Moderators (coefficients 1:24):
## QM(df = 24) = 147.1773, p-val < .0001
##
## Model Results:
##
## estimate
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment 0.3076
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 0.0876
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment 0.4751
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment -0.0980
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.0299
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 0.1172
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.7433
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1154
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution -1.0873
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution -0.2194
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 0.0173
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution -0.0897
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.1289
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 0.2036
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -0.2680
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -0.2483
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource -0.1147
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource -0.8293
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource -0.5378
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource -0.1989
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -0.7288
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -0.7240
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.0670
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource -0.2270
## se
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment 0.1629
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 0.2250
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment 0.1009
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment 0.2509
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.3119
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 0.5969
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1019
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.1986
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution 0.4160
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution 0.3152
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 0.1640
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution 0.1999
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.3957
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 0.3577
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.1771
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 0.1963
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource 0.2557
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource 0.9515
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource 0.2155
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource 0.3399
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource 0.2431
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 0.3619
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.2220
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource 0.4022
## zval
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment 1.8881
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 0.3895
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment 4.7067
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment -0.3907
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.0959
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 0.1964
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -7.2947
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.5809
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution -2.6135
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution -0.6962
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 0.1053
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution -0.4488
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.3259
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 0.5691
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -1.5130
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -1.2647
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource -0.4487
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource -0.8716
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource -2.4961
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource -0.5853
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -2.9984
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -2.0003
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.3019
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource -0.5644
## pval
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment 0.0590
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 0.6969
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment <.0001
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment 0.6960
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.9236
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 0.8443
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment <.0001
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.5613
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution 0.0090
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution 0.4863
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 0.9161
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution 0.6536
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.7445
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 0.5693
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.1303
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 0.2060
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource 0.6537
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource 0.3834
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource 0.0126
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource 0.5583
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource 0.0027
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 0.0455
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.7627
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource 0.5725
## ci.lb
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment -0.0117
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment -0.3534
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment 0.2773
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment -0.5898
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment -0.5815
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment -1.0526
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.9430
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.2739
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution -1.9027
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution -0.8371
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution -0.3041
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution -0.4814
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution -0.9045
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution -0.4975
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution -0.6151
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution -0.6331
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource -0.6159
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource -2.6941
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource -0.9601
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource -0.8651
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -1.2052
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -1.4334
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource -0.3680
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource -1.0152
## ci.ub
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment 0.6268
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 0.5286
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment 0.6730
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment 0.3937
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 0.6412
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 1.2871
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment -0.5436
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 0.5047
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution -0.2719
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution 0.3983
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 0.3386
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution 0.3020
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 0.6466
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 0.9046
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 0.0792
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 0.1365
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource 0.3864
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource 1.0356
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource -0.1155
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource 0.4672
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource -0.2524
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource -0.0146
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 0.5020
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource 0.5613
##
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment .
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment ***
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment ***
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution **
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource *
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource **
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource *
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest.default(
x = Q2.R1$beta,
sei = Q2.R1$se,
ci.lb = Q2.R1$ci.lb,
ci.ub = Q2.R1$ci.ub,
annotate = TRUE,
showweights = T,
header = F,
slab = c(
"AQ:EE:prevalence",
"TE:EE:prevalence",
"AQ:EE:intensity",
"TE:EE:intensity",
"AQ:EE:fecundity",
"TE:EE:fecundity",
"AQ:EE:survivorship",
"TE:EE:survivorship",
"AQ:CP:prevalence",
"TE:CP:prevalence",
"AQ:CP:intensity",
"TE:CP:intensity",
"AQ:CP:fecundity",
"TE:CP:fecundity",
"AQ:CP:survivorship",
"TE:CP:survivorship",
"AQ:RL:prevalence",
"TE:RL:prevalence",
"AQ:RL:intensity",
"TE:RL:intensity",
"AQ:RL:fecundity",
"TE:RL:fecundity",
"AQ:RL:survivorship",
"TE:RL:survivorship"
)
)
par_short <- c("AQ:EE:prevalence",
"TE:EE:prevalence",
"AQ:EE:intensity",
"TE:EE:intensity",
"AQ:EE:fecundity",
"TE:EE:fecundity",
"AQ:EE:survivorship",
"TE:EE:survivorship",
"AQ:CP:prevalence",
"TE:CP:prevalence",
"AQ:CP:intensity",
"TE:CP:intensity",
"AQ:CP:fecundity",
"TE:CP:fecundity",
"AQ:CP:survivorship",
"TE:CP:survivorship",
"AQ:RL:prevalence",
"TE:RL:prevalence",
"AQ:RL:intensity",
"TE:RL:intensity",
"AQ:RL:fecundity",
"TE:RL:fecundity",
"AQ:RL:survivorship",
"TE:RL:survivorship")
pal <- RColorBrewer::brewer.pal(11,"Spectral")
Q2.R1_OP <-
orchard_plot(
Q2.R1,
xlab = "Standardised mean difference (Hedge's g)",
transfm = "none",
angle = 0,
alpha = 0.7
) +
theme(legend.position = c(0.3, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
scale_fill_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 4)) +
scale_colour_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 4)) +
scale_y_discrete(labels = par_short)## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.Q2.R1_OP
Make Table with results from Q2.R1
TableQ2.R1 <-
data.frame(
Environment = rep(c(
"Aquatic",
"Terrestrial"), 12),
Gradient.category = rep(
c(
"Environment",
"Pollution",
"Resource"
),
each = 8
),
Trait.type = rep(c(
"Prevalence",
"Prevalence",
"Intensity",
"Intensity",
"Fecundity",
"Fecundity",
"Survivorship",
"Survivorship"), 3
),
Overall_mean = round(Q2.R1$beta, 3),
Lower_95 = round(Q2.R1$ci.lb, 3),
Upper_95 = round(Q2.R1$ci.ub, 3),
P_value = round(Q2.R1$pval, 3))
N_effects <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_effects = length(ID))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.N_experiments <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_experiments = length(unique(Experiment)))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.N_hosts <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_hosts = length(unique(Host)))## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.TableQ2.R1 <-
left_join(TableQ2.R1, N_effects) %>% left_join(N_experiments) %>% left_join(N_experiments) %>%
dplyr::rename(Stressor_type = Gradient.category,
Response_trait = Trait.type) %>%
mutate(Stressor_type = plyr::revalue(
Stressor_type,
c(
"Environment" = "Endogenous Environment",
"Pollution" = "Chemical Pollution",
"Resource" = "Resource Limitation"
)
))## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type,
## N_experiments)`rownames(TableQ2.R1) <- NULL
TableQ2.R1 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Environment | Stressor_type | Response_trait | Overall_mean | Lower_95 | Upper_95 | P_value | N_effects | N_experiments |
|---|---|---|---|---|---|---|---|---|
| Aquatic | Endogenous Environment | Prevalence | 0.308 | -0.012 | 0.627 | 0.059 | 41 | 21 |
| Terrestrial | Endogenous Environment | Prevalence | 0.088 | -0.353 | 0.529 | 0.697 | 35 | 7 |
| Aquatic | Endogenous Environment | Intensity | 0.475 | 0.277 | 0.673 | 0.000 | 135 | 50 |
| Terrestrial | Endogenous Environment | Intensity | -0.098 | -0.590 | 0.394 | 0.696 | 17 | 5 |
| Aquatic | Endogenous Environment | Fecundity | 0.030 | -0.581 | 0.641 | 0.924 | 7 | 3 |
| Terrestrial | Endogenous Environment | Fecundity | 0.117 | -1.053 | 1.287 | 0.844 | 2 | 1 |
| Aquatic | Endogenous Environment | Survivorship | -0.743 | -0.943 | -0.544 | 0.000 | 125 | 56 |
| Terrestrial | Endogenous Environment | Survivorship | 0.115 | -0.274 | 0.505 | 0.561 | 37 | 10 |
| Aquatic | Chemical Pollution | Prevalence | -1.087 | -1.903 | -0.272 | 0.009 | 8 | 4 |
| Terrestrial | Chemical Pollution | Prevalence | -0.219 | -0.837 | 0.398 | 0.486 | 11 | 5 |
| Aquatic | Chemical Pollution | Intensity | 0.017 | -0.304 | 0.339 | 0.916 | 48 | 17 |
| Terrestrial | Chemical Pollution | Intensity | -0.090 | -0.481 | 0.302 | 0.654 | 26 | 13 |
| Aquatic | Chemical Pollution | Fecundity | -0.129 | -0.905 | 0.647 | 0.745 | 7 | 4 |
| Terrestrial | Chemical Pollution | Fecundity | 0.204 | -0.498 | 0.905 | 0.569 | 9 | 3 |
| Aquatic | Chemical Pollution | Survivorship | -0.268 | -0.615 | 0.079 | 0.130 | 46 | 18 |
| Terrestrial | Chemical Pollution | Survivorship | -0.248 | -0.633 | 0.136 | 0.206 | 24 | 14 |
| Aquatic | Resource Limitation | Prevalence | -0.115 | -0.616 | 0.386 | 0.654 | 14 | 4 |
| Terrestrial | Resource Limitation | Prevalence | -0.829 | -2.694 | 1.036 | 0.383 | 1 | 1 |
| Aquatic | Resource Limitation | Intensity | -0.538 | -0.960 | -0.116 | 0.013 | 25 | 9 |
| Terrestrial | Resource Limitation | Intensity | -0.199 | -0.865 | 0.467 | 0.558 | 10 | 4 |
| Aquatic | Resource Limitation | Fecundity | -0.729 | -1.205 | -0.252 | 0.003 | 18 | 6 |
| Terrestrial | Resource Limitation | Fecundity | -0.724 | -1.433 | -0.015 | 0.045 | 13 | 2 |
| Aquatic | Resource Limitation | Survivorship | 0.067 | -0.368 | 0.502 | 0.763 | 21 | 7 |
| Terrestrial | Resource Limitation | Survivorship | -0.227 | -1.015 | 0.561 | 0.573 | 6 | 3 |
Does accounting for the transmission environment reduce heterogeneity? We calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix.
Q2.R1.R <-
rma.mv(
g ~ Environment:Trait.type:Gradient.category -1 ,
V = dat_2$var.g,
random = list( ~ 1 | Experiment, ~ 1 | ID),
data = dat_2,
method = "REML"
)
I2.Q2.R1 <- I2(Q2.R1.R, v = dat_1$var.g, obs = "ID")*100## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'I2.Q2.R1## I2_Est. 2.5% CI 97.5% CI
## Experiment 29.89 24.21 35.83
## total 94.35 93.74 94.89
Comment 1.R1 - Do infected and uninfected hosts differ in fitness?
We’ll address this question by Reviewer 1 in the control condition (i.e. in the absence of stressors)
Re arrange data
Calculate OR and approximate variance from contingency tables
Calculate SD for comparisons with continuous normally distributed variables
Get effect sizes and variances